title: non linear regression with random effect and lsoda
body: <p>I am facing a problem I do not manage to solve. I would like to use <code>nlme</code> or <code>nlmODE</code> to perform a non linear regression with random effect using as a model the solution of a second order differential equation with fixed coefficients (a damped oscillator).</p>
  <p>I manage to use <code>nlme</code> with simple models, but it seems that the use of <code>deSolve</code> to generate the solution of the differential equation causes a problem. Below an example, and the problems I face.</p>
  <h1>The data and functions</h1>
  <p>Here is the function to generate the solution of the differential equation using <code>deSolve</code>:</p>
  <pre><code>library(deSolve)
  ODE2_nls &lt;- function(t, y, parms) {
    S1 &lt;- y[1]
    dS1 &lt;- y[2]
    dS2 &lt;- dS1
    dS1 &lt;- - parms[&quot;esp2omega&quot;]*dS1  - parms[&quot;omega2&quot;]*S1 + parms[&quot;omega2&quot;]*parms[&quot;yeq&quot;]
    res &lt;- c(dS2,dS1)
    list(res)}
  
  solution_analy_ODE2 = function(omega2,esp2omega,time,y0,v0,yeq){
    parms  &lt;- c(esp2omega = esp2omega,
                omega2 = omega2,
                yeq = yeq)
    xstart = c(S1 =  y0, dS1 = v0)
    out &lt;-  lsoda(xstart, time, ODE2_nls, parms)
    return(out[,2])
  }
  </code></pre>
  <p>I can generate a solution for a given period and damping factor, as for example here a period of 20 and a slight damping of 0.2:</p>
  <pre class="lang-r prettyprint-override"><code>
  # small example:
  time &lt;- 1:100
  period &lt;- 20 # period of oscillation
  amort_factor &lt;- 0.2
  omega &lt;- 2*pi/period # agular frequency
  oscil &lt;- solution_analy_ODE2(omega^2,amort_factor*2*omega,time,1,0,0)
  plot(time,oscil)
  
  </code></pre>
  <p><a href="https://i.stack.imgur.com/vEwIO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/vEwIO.png" alt="enter image description here" /></a></p>
  <p>Now I generate a panel of 10 individuals with a random starting phase (i.e. different starting position and velocity). The goal is to perform a non linear regression with <strong>random effect on the starting values</strong></p>
  <pre><code>library(data.table)
  # generate panel
  Npoint &lt;- 100 # number of time poitns
  Nindiv &lt;- 10 # number of individuals
  period &lt;- 20 # period of oscillation
  amort_factor &lt;- 0.2
  omega &lt;- 2*pi/period # agular frequency
  # random phase
  phase &lt;- sample(seq(0,2*pi,0.01),Nindiv)
  # simu data:
  data_simu &lt;- data.table(time = rep(1:Npoint,Nindiv), ID = rep(1:Nindiv,each = Npoint))
  
  # signal generation
  data_simu[,signal := solution_analy_ODE2(omega2 = omega^2,
                                           esp2omega = 2*0.2*omega,
                                           time = time,
                                           y0 = sin(phase[.GRP]),
                                           v0 = omega*cos(phase[.GRP]),
                                           yeq = 0)+ 
              rnorm(.N,0,0.02),by = ID]
  </code></pre>
  <p>If we have a look, we have a proper dataset:</p>
  <pre><code>library(ggplot2)
  ggplot(data_simu,aes(time,signal,color = ID))+
    geom_line()+
    facet_wrap(~ID)
  </code></pre>
  <p><a href="https://i.stack.imgur.com/1EFtO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/1EFtO.png" alt="enter image description here" /></a></p>
  <h1>The problems</h1>
  <h2>Using nlme</h2>
  <p>Using <code>nlme</code> with similar syntax working on simpler examples (non linear functions not using deSolve), I tried:</p>
  <pre><code>fit &lt;- nlme(model = signal ~ solution_analy_ODE2(esp2omega,omega2,time,y0,v0,yeq), 
       data = data_simu,
       fixed = esp2omega + omega2 + y0 + v0 + yeq ~ 1,
       random = y0 ~ 1 ,
       groups = ~ ID, 
       start = c(esp2omega = 0.08, 
                 omega2 = 0.04,
                 yeq = 0,
                 y0 = 1,
                 v0 = 0))
  </code></pre>
  <p>I obtain:</p>
  <blockquote>
  <p>Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (2) must equal the length of the initial conditions vector (2000)</p>
  </blockquote>
  <p>The traceback:</p>
  <pre><code>12. stop(paste(&quot;The number of derivatives returned by func() (&quot;, length(tmp[[1]]), &quot;) must equal the length of the initial conditions vector (&quot;, length(y), &quot;)&quot;, sep = &quot;&quot;))
  11. checkFunc(Func2, times, y, rho)
  10. lsoda(xstart, time, ODE2_nls, parms)
  9. solution_analy_ODE2(omega2, esp2omega, time, y0, v0, yeq)
  .
  .
  </code></pre>
  <p>I looks like <code>nlme</code> is trying to pass a vector of starting condition to <code>solution_analy_ODE2</code>, and causes an error in <code>checkFunc</code> from <code>lasoda</code>.</p>
  <p>I tried using <code>nlsList</code>:</p>
  <pre><code>test &lt;- nlsList(model = signal ~ solution_analy_ODE2(omega2,esp2omega,time,y0,v0,yeq) | ID, 
          data = data_simu, 
          start = list(esp2omega = 0.08, omega2 = 0.04,yeq = 0,
                       y0 = 1,v0 = 0),
          control = list(maxiter=150, warnOnly=T,minFactor = 1e-10), 
          na.action = na.fail, pool = TRUE)
  head(test)
  </code></pre>
  <pre><code>Call:
    Model: signal ~ solution_analy_ODE2(omega2, esp2omega, time, y0, v0, yeq) | ID 
     Data: data_simu 
  
  Coefficients:
     esp2omega     omega2           yeq         y0          v0
  1  0.1190764 0.09696076  0.0007577956 -0.1049423  0.30234654
  2  0.1238936 0.09827158 -0.0003463023  0.9837386  0.04773775
  3  0.1280399 0.09853310 -0.0004908579  0.6051663  0.25216134
  4  0.1254053 0.09917855  0.0001922963 -0.5484005 -0.25972829
  5  0.1249473 0.09884761  0.0017730823  0.7041049  0.22066652
  6  0.1275408 0.09966155 -0.0017522320  0.8349450  0.17596648
  </code></pre>
  <p>We can see that te non linear fit works well on individual signals. Now if I want to perform a regression of the dataset with random effects, the syntax should be:</p>
  <pre><code>fit &lt;- nlme(test, 
       random = y0 ~ 1 ,
       groups = ~ ID, 
       start = c(esp2omega = 0.08, 
                 omega2 = 0.04,
                 yeq = 0,
                 y0 = 1,
                 v0 = 0))
  </code></pre>
  <p>But I obtain the exact same error message.</p>
  <p>I then tried using <code>nlmODE</code>, following Bne Bolker's comment on <a href="https://stackoverflow.com/questions/56254840/fitting-first-order-equation-with-nlme-and-lsoda">a similar question I asked some years ago</a></p>
  <h2>using nlmODE</h2>
  <pre><code>library(nlmeODE)
  datas_grouped &lt;- groupedData( signal ~ time | ID, data = data_simu, 
                                labels = list (x = &quot;time&quot;, y = &quot;signal&quot;), 
                                units = list(x =&quot;arbitrary&quot;, y = &quot;arbitrary&quot;))
  
  modelODE &lt;- list( DiffEq = list(dS2dt = ~ S1,
                                  dS1dt = ~ -esp2omega*S1  - omega2*S2 + omega2*yeq),
                    ObsEq = list(yc = ~ S2),
                    States = c(&quot;S1&quot;,&quot;S2&quot;),
                    Parms = c(&quot;esp2omega&quot;,&quot;omega2&quot;,&quot;yeq&quot;,&quot;ID&quot;), 
                    Init = c(y0 = 0,v0 = 0))
  
  resnlmeode = nlmeODE(modelODE, datas_grouped) 
  assign(&quot;resnlmeode&quot;, resnlmeode, envir = .GlobalEnv)
  #Fitting with nlme the resulting function
  model &lt;- nlme(signal ~ resnlmeode(esp2omega,omega2,yeq,time,ID), 
                data = datas_grouped, 
                fixed = esp2omega + omega2 + yeq + y0 + v0  ~ 1, 
                random = y0 + v0 ~1,
                start = c(esp2omega = 0.08, 
                          omega2 = 0.04,
                          yeq = 0,
                          y0 = 0,
                          v0 = 0)) # 
  
  </code></pre>
  <p>I get the error:</p>
  <blockquote>
  <p>Error in resnlmeode(esp2omega, omega2, yeq, time, ID) : object 'yhat' not found</p>
  </blockquote>
  <p>Here I don't understand where the error comes from, nor how to solve it.</p>
  <h1>Questions</h1>
  <ul>
  <li>Can you reproduce the problem ?</li>
  <li>Does anyone have an idea to solve this problem, using either <code>nlme</code> or <code>nlmODE</code> ?</li>
  <li>If not, is there a solution using an other package ? I saw <code>nlmixr</code> (<a href="https://cran.r-project.org/web/packages/nlmixr/index.html" rel="nofollow noreferrer">https://cran.r-project.org/web/packages/nlmixr/index.html</a>), but I don't know it, the instalation is complicated and it was recently remove from CRAN</li>
  </ul>
  <hr />
  <h1>Edits</h1>
  <p>@tpetzoldt suggested a nice way to debug <code>nlme</code> behavior, and it surprised me a lot. Here is a working example with a non linear function, where I generate a set of 5 individual with a random parameter varying between individuals :</p>
  <pre><code>reg_fun = function(time,b,A,y0){
    cat(&quot;time : &quot;,length(time),&quot; b :&quot;,length(b),&quot; A : &quot;,length(A),&quot; y0: &quot;,length(y0),&quot;\n&quot;)
    out &lt;- A*exp(-b*time)+(y0-1)
    cat(&quot;out : &quot;,length(out),&quot;\n&quot;)
    tmp &lt;- cbind(b,A,y0,time,out)
    cat(apply(tmp,1,function(x) paste(paste(x,collapse = &quot; &quot;),&quot;\n&quot;)),&quot;\n&quot;)
    return(out)
  }
  
  time &lt;- 0:10*10
  ramdom_y0 &lt;- sample(seq(0,1,0.01),10)
  Nid &lt;- 5
  data_simu &lt;- 
  data.table(time = rep(time,Nid),
             ID = rep(LETTERS[1:Nid],each = length(time)) )[,signal := reg_fun(time,0.02,2,ramdom_y0[.GRP]) + rnorm(.N,0,0.1),by = ID]
  </code></pre>
  <p>The cats in the function give here:</p>
  <pre><code>time :  11  b : 1  A :  1  y0:  1 
  out :  11 
  0.02 2 0.64 0 1.64 
   0.02 2 0.64 10 1.27746150615596 
   0.02 2 0.64 20 0.980640092071279 
   0.02 2 0.64 30 0.737623272188053 
   0.02 2 0.64 40 0.538657928234443 
   0.02 2 0.64 50 0.375758882342885 
   0.02 2 0.64 60 0.242388423824404 
   0.02 2 0.64 70 0.133193927883213 
   0.02 2 0.64 80 0.0437930359893108 
   0.02 2 0.64 90 -0.0294022235568269 
   0.02 2 0.64 100 -0.0893294335267746
  .
  .
  .
  </code></pre>
  <p>Now I do with <code>nlme</code>:</p>
  <pre><code>nlme(model = signal ~ reg_fun(time,b,A,y0), 
       data = data_simu,
       fixed = b + A + y0 ~ 1,
       random = y0 ~ 1 ,
       groups = ~ ID, 
       start = c(b = 0.03, A = 1,y0 = 0))
  </code></pre>
  <p>I get:</p>
  <pre><code>time :  55  b : 55  A :  55  y0:  55 
  out :  55 
  0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   
  time :  55  b : 55  A :  55  y0:  55 
  out :  55 
  0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
   0.03 1 0 0 0 
   0.03 1 0 10 -0.259181779318282 
   0.03 1 0 20 -0.451188363905974 
   0.03 1 0 30 -0.593430340259401 
   0.03 1 0 40 -0.698805788087798 
   0.03 1 0 50 -0.77686983985157 
   0.03 1 0 60 -0.834701111778413 
   0.03 1 0 70 -0.877543571747018 
   0.03 1 0 80 -0.909282046710588 
   0.03 1 0 90 -0.93279448726025 
   0.03 1 0 100 -0.950212931632136 
  ...
  </code></pre>
  <p>So <code>nlme</code> binds 5 time (the number of individual) the time vector and pass it to the function, with the parameters repeated the same number of time. Which is of course not compatible with the way <code>lsoda</code> and my function works.</p>