[R-sig-ME] nlme: "Singularity in backsolve" error in ODE-defined, but not in closed-form model

[Tobias Sing]

Dear all,

I would appreciate very much some help/advice on a problem related to
using nlme together with odesolve. I am aware that the issue might
also be related to odesolve rather than to nlme. Still, the error
message comes from nmle, and I know several people working with nlme
who have run into similar problems, without knowing a solution
(platform: nlme 3.1-79, R 2.4.1, Win XP).

I have been trying to come up with a minimalistic example of using
analytical versus ODE-defined models in nlme. Specifically, I am
sampling artificial data from the model
y_ij = exp(-(beta + b_i) * t) + eps_ij
(i = number of individuals, j=observations for an individual),
which corresponds to the ODE dy/dt = -beta * y(t).

As in the examples in "the book", I am fitting three models:
1.) Pooling the data from all individuals
2.) Separately for all individuals
3.) Mixed-effects model

All three steps work fine for the _analytically_ defined models (and
the estimates are all good), but when I use the equivalent
_ODE_-defined models, only the first two steps are fine, and the call
to nlme (for fitting model 3) always stops with

--> "Error in MEEM(...): Singularity in backsolve at level 0, block 1)".

Here is what I tried to solve this problem, without success:
a.) Setting atol and rtol in odesolve as low as possible to be
virtually identical to the output of the analytical expression. In R's
output, I don't see any remaining differences. For example, I have
analytical(TIME=10, beta=0.75) = ode(TIME=10, beta=0.75) =
0.0005530844
b.) Playing around with the nlme control options msVerbose, tolerance,
msTol, pnlsTol, and maxiter.
c.) Googling for possible solutions

I would appreciate any help/advice to solve this problem.
Many thanks,
  Tobias

P.S. Here is the full code to reproduce the problem:

## sample artificial exponential decay data with noise:
library(nlme)
library(odesolve)
beta <- 1
sigma <- 0.025
sigma.b <- 0.2
times <- seq(0,10, by=0.5)

n.individuals <- 10
b <- rnorm(n.individuals, 0, sigma.b)
eps <- matrix(rnorm(n.individuals * length(times), 0, sigma),
length(times), n.individuals)

my.data <- as.list(1:n.individuals)
for (i in 1:n.individuals) {
  my.data[[i]] <- data.frame(ID=rep(i, length(times)),
                       TIME=times,
                       DV=exp(- (beta + b[i]) * times) + eps[, i])
  my.data[[i]][ my.data[[i]] <= 0] <- 10^(-5)
}
my.data <- do.call(rbind, my.data)

## fit analytical models (pooled, individual, and mixed):
## EVERYTHING WORKS FINE HERE!
analytical <- function(TIME, beta) {
  exp(-beta * TIME)
}
nls(DV ~ analytical(TIME, beta), data=my.data, start=list(beta=2))
model.list <- nlsList(DV ~ analytical(TIME, beta) | ID, data=my.data,
start=list(beta=2))
model.nlme <- nlme(model.list)

## fit equivalent ODE-defined models (pooled, individual, and mixed):
## --> ERROR "Singularity in backsolve" for mixed model
yprime <- function(t, y, parms) {
  dydt <- -parms['beta'] * y
  list(dydt)
}
ode <- function(TIME, beta) {
  y0 <- 1
  ans <- lsoda(y0, c(0,TIME), yprime, c(beta=beta), atol=10^(-10),
rtol=10^(-10))
  ans[2,2]
}
ode.vectorized <- Vectorize(ode, 'TIME')

nls(DV ~ ode.vectorized(TIME, beta), data=my.data, start=list(beta=2))
model.list <- nlsList(DV ~ ode.vectorized(TIME, beta) | ID,
data=my.data, start=list(beta=2))
model.nlme <- nlme(model.list)

"Error in MEEM(object, conLin, control$niterEM) :
        Singularity in backsolve at level 0, block 1"