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))

Rathe

have you tried the nlmeODE package? – Aspa

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 – Rathe

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 ... – Aspa

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. – Ethicize

@Ethicize 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 – Rathe

@Rathe this will probably be helpful: #40025639 – Ethicize

Answer 1

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

Aspa

I made with two subjects as an example. The goal is to fit to a real panel of (simulated) data, e.g. 100 individuals. – Rathe

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. – Rathe