fitting first order equation with nlme and lsoda

Question

I a trying to fit a first order differential model using nlme and lsoda. Here is the basic idea: I first define the function allowing to generate the solution of the differential equation:

library(deSolve)

ODE1 <- function(time, x, parms) {with(as.list(c(parms, x)), {
  import <- excfunc(time)
  dS <- import*k/tau - (S-yo)/tau 
  res <- c(dS)
  list(res)})}


solution_ODE1 = function(tau1,k1,yo1,excitation,time){
  excfunc <- approxfun(time, excitation, rule = 2)
  parms  <- c(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
  xstart = c(S = yo1)
  out <-  lsoda(xstart, time, ODE1, parms)
  return(out[,2])
}

I then generate data following the equation for two IDs:

time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(signal = c(solution_ODE1(3,2,0.1,excitation,time)+rnorm(length(time),0,0.1),
                                   solution_ODE1(3.2,1.5,0.3,excitation,time)+rnorm(length(time),0,0.1)),
                        time = rep(time,2),
                        excitation = rep(excitation,2),
                        ID = rep(c("A","B"),each = length(time)))

Here is what it looks like :

library(ggplot2)
ggplot(simu_data)+
  geom_point(aes(time,signal,color = "signal"),size = 2)+
  geom_line(aes(time,excitation,color = "excitation"))+
  facet_wrap(~ID)

I am then trying to fit using nlme:

fit1 <- nlme(signal ~ solution_ODE1(damping,gain,eq,excitation,time),
             data = simu_data,
             fixed = damping + gain + eq ~1,
             random =  damping   ~ 1 ,
             groups = ~ ID,
             start = c(damping = 5, gain = 1,eq = 0))

I am getting this eror, that I don't get:

Error in eval(substitute(expr), data, enclos = parent.frame()) : object 'k' not found

The traceback shows that the error comes from the ODE1 model, which works when generating values.

16.    eval(substitute(expr), data, enclos = parent.frame()) 
15.    eval(substitute(expr), data, enclos = parent.frame()) 
14.    with.default(as.list(c(parms, x)), {
    import <- excfunc(time)
    dS <- import * k/tau - (S - yo)/tau
    res <- c(dS) ... 
13.    with(as.list(c(parms, x)), {
    import <- excfunc(time)
    dS <- import * k/tau - (S - yo)/tau
    res <- c(dS) ... 
12.    func(time, state, parms, ...) 
11.    Func2(times[1], y) 
10.    eval(Func2(times[1], y), rho) 
9.    checkFunc(Func2, times, y, rho) 
8.    lsoda(xstart, time, ODE1, parms) 
7.    solution_ODE1(damping, gain, eq, excitation, time) 
6.    eval(model, data.frame(data, pars)) 
5.    eval(model, data.frame(data, pars)) 
4.    eval(modelExpression[[2]], envir = nlEnv) 
3.    eval(modelExpression[[2]], envir = nlEnv) 
2.    nlme.formula(signal ~ solution_ODE1(damping, gain, eq, excitation, 
    time), data = simu_data, fixed = damping + gain + eq ~ 1, 
    random = damping ~ 1, groups = ~ID, start = c(damping = 5, 
        gain = 1, eq = 0)) 
1.    nlme(signal ~ solution_ODE1(damping, gain, eq, excitation, time), 
    data = simu_data, fixed = damping + gain + eq ~ 1, random = damping ~ 
        1, groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))

Does anyone have an idea How I should proceed ?

Edit

I tried to modify following the advise of mikeck:

ODE1 <- function(time, x, parms) {
  import <- parms$excfunc(time)
  dS <- import*parms$k/parms$tau - (x["S"]-parms$yo)/parms$tau 
  res <- c(dS)
  list(res)}

Generating the data works without problems. But use of nlme gives now:

Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (0) must equal the length of the initial conditions vector (100)

with the following traceback:

> traceback()
10: stop(paste("The number of derivatives returned by func() (", 
        length(tmp[[1]]), ") must equal the length of the initial conditions vector (", 
        length(y), ")", sep = ""))
9: checkFunc(Func2, times, y, rho)
8: lsoda(xstart, time, ODE1, parms) at #5
7: solution_ODE1(damping, gain, eq, excitation, time)
6: eval(model, data.frame(data, pars))
5: eval(model, data.frame(data, pars))
4: eval(modelExpression[[2]], envir = nlEnv)
3: eval(modelExpression[[2]], envir = nlEnv)
2: nlme.formula(signal ~ solution_ODE1(damping, gain, eq, excitation, 
       time), data = simu_data, fixed = damping + gain + eq ~ 1, 
       random = damping ~ 1, groups = ~ID, start = c(damping = 5, 
           gain = 1, eq = 0))
1: nlme(signal ~ solution_ODE1(damping, gain, eq, excitation, time), 
       data = simu_data, fixed = damping + gain + eq ~ 1, random = damping ~ 
           1, groups = ~ID, start = c(damping = 5, gain = 1, eq = 0))

I am actually trying. I am having a bit of hard time with it, but maybe it will do the trick. I am still happy to get a solution/explanation for this weird behavior — denis, May 27 '19 at 15:07
I've made a few adjustments - parms should be use list(), not c(), and I've made `xstart <- yo1` (then referring directly to `x` in `ODE1`, but I'm stilling getting an "illegal input" message ... — Ben Bolker, May 28 '19 at 01:46
Have you tried redefining `ODE1()` to not use `with()`, i.e. instead use `parms$k` etc.? The error message looks like it might be a scoping issue that is cropping up somehow. — mikeck, May 30 '19 at 23:36
@mikeck I tried, it changed the error message. I edited my question. I don't get what `nlme` does internally, but it looks it give vector of initial condition to the function, creating thus errors — denis, Jun 03 '19 at 08:52
@denis this will probably be helpful: https://stackoverflow.com/questions/40025139/solving-ode-with-desolve-in-r-number-of-derivatives-error — mikeck, Jun 03 '19 at 16:39

score 3 · Answer 1 · answered May 28 '19 at 01:51

In your example, your times vector doesn't run monotonically. I think that messes with lsoda. What is the context/meaning of the way that time works here? It doesn't really make sense to fit a random-effects model with two groups. Are you trying to fit the same curve to two independent time series?

Here's a stripped-down example, with some adjustments (not everything can be collapsed to a numeric vector without losing necessary structure):

library(deSolve)
ODE1 <- function(time, x, parms) {
    with(as.list(parms), {
        import <- excfunc(time)
        dS <- import*k/tau - (x-yo)/tau 
        res <- c(dS)
        list(res)
    })
}
solution_ODE1 = function(tau1,k1,yo1,excitation,time){
    excfunc <- approxfun(time, excitation, rule = 2)
    parms  <- list(tau = tau1, k = k1, yo = yo1, excfunc = excfunc)
    xstart = yo1
    out <-  lsoda(xstart, time, ODE1, parms)
    return(out[,2])
}
time <- 0:49
excitation <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10),rep(0,10))
simu_data <- data.frame(time = rep(time,2),
                        excitation = rep(excitation,2))
svec <- c(damping = 3, gain = 1.75, eq = 0.2)

This works:

with(c(simu_data, as.list(svec)),
     solution_ODE1(damping,gain,eq,excitation[1:50],time[1:50]))

But if we include one more step (so that time resets to 0), it fails:

with(c(simu_data, as.list(svec)),
     solution_ODE1(damping,gain,eq,excitation[1:51],time[1:51]))

Error in lsoda(xstart, time, ODE1, parms) : illegal input detected before taking any integration steps - see written message

I made with two subjects as an example. The goal is to fit to a real panel of (simulated) data, e.g. 100 individuals. — denis, May 28 '19 at 08:05
The time vector run monotically for each subject (`ID` variable in `simu_data`). Of course if you include time values from the second subject into lsoda, it will cause an error, because you will have twice the same time value. Here lsoda works perfectly to generate the data, as I use it to produce the `simu_data` dataset. But it seems that `nlme` do something while varying the parameters that cause an error, and I can't figure out what. — denis, May 28 '19 at 08:10

fitting first order equation with nlme and lsoda

Edit

1 Answers1

Linked