"Collection operator c error" when calling OpenBUGS from R

Question

I'm trying to use R2OpenBUGS to run a differential equation solver in OpenBUGS. I have tried this with the exponential decay example (Example01.odc) in the OpenBUGS/Diff/Examples folder.

The first part of this example is a simulation, the second is the inference version, which I try to apply. Both run OK directly within OpenBUGS.

When I try to run it from R, I get an error in the log :

expected the collection operator c error pos 1443
variable ngrid is not defined 

I think the second is a follow up from the first.

Below is what I have coded to run it from R :

# example of adjusting a differential equation using OpenBUGS   
library(R2OpenBUGS)
library(coda)
workingdir <- "M:/data/R/OpenBUGSDiff"
setwd(workingdir)
# time points
tgrid <- c(0.00000, 0.20000, 0.40000, 0.60000, 0.80000,
       1.00000, 1.20000, 1.40000, 1.60000, 1.80000,
       2.00000, 2.20000, 2.40000, 2.60000, 2.80000,
       3.00000, 3.20000, 3.40000, 3.60000, 3.80000,
       4.00000, 4.20000, 4.40000, 4.60000, 4.80000,
       5.00000, 5.20000, 5.40000, 5.60000, 5.80000,
       6.00000, 6.20000, 6.40000, 6.60000, 6.80000,
       7.00000, 7.20000, 7.40000, 7.60000, 7.80000,
       8.00000, 8.20000, 8.40000, 8.60000, 8.80000,
       9.00000, 9.20000, 9.40000, 9.60000, 9.80000,
       10.00000)
obs <- c(  0.7887,1.241,0.7051,0.7388,0.3903,
    0.2632,0.1174,0.549,-0.1767,0.02938,
    0.154,0.1465,0.07878,-0.5569,0.01293,
    0.2905,-0.2665,-0.3881,0.02517,-0.138,
    0.4004,0.2859,-0.1217,0.3961,0.3813,
    0.1846,-0.3581,0.3293,0.04089,0.01972,
    0.3203,0.5294,-0.1389,-0.3732,0.1341,
    -0.02432,0.2261,-0.3612,0.3131,-0.258,
    0.02948,-0.0208,0.1066,0.3796,-0.2645,
    0.1035,0.1001,-0.2415,0.06739,-0.1554,
    -0.2388)
# inference model 
Modele <- function() 
{
  solution[1:ngrid, 1] <- 
    ode(init[1], 
        tgrid[1:ngrid], 
        D(C[1], t), 
        origin, 
        tol)
  D(C[1], t) <- -lambda * C[1] 
  log.lambda ~ dnorm(0.0, tau.lambda); 
  lambda <- exp(log.lambda)
  for (i in 1:ngrid) 
  {
    obs[i] ~ dnorm(solution[i, 1], tau) 
  }
  tau ~ dgamma(a, b)
}
write.model(Modele,"Diff.txt")
# data definition
origin = 0.0
tol = 1.0E-3
ngrid = 51 
init = c(1.0)
a = 0.001
b = 0.001
tau.lambda = 0.01
data  <- list(obs=obs,
          tgrid=tgrid,
          origin=origin,
          tol=tol,
          ngrid=ngrid,
          init=init,
          a=a,
          b=b,
          tau.lambda=tau.lambda)
inits <- function(){
  list(log.lambda = 1.0,tau = 0.01)
}
diff.inf <- bugs(data,inits,model.file = "Diff.txt",
              parameters = c("lambda","tau"),
              n.chains = 1, n.iter = 3000,n.burnin=500,n.thin=10,
              working.directory=workingdir,
              debug = TRUE,
              codaPkg=T)

And here is the complete error message from BUGS :

model is syntactically correct
expected the collection operator c error pos 1443
variable ngrid is not defined
model must have been compiled but not updated to be able to change RN generator
BugsCmds:NoCompileInits
model must be compiled before generating initial values
model must be initialized before updating
model must be initialized before monitors used
model must be initialized before monitors used
model must be initialized before monitors used
model must be initialized before DIC can be monitored
model must be initialized before updating
model must be initialized before monitors used
DIC monitor not set

I have tried to replace the "C[]" in the model by "obs" as I don't see where the "C" is defined, but this did not change anything.

I checked for "e" instead of "E" in the scientific number notation, which has been suggested elsewhere to create "variable not defined " errors in BUGS, but that does not seem to be the problem.

Any ideas are welcome, or any other code that shows how to use "solution" or "ode" from within R.

Thanks in advance

Oliver

Answer 1

There is a BUGS error related to the init node. In your model above you have an index associated with it (init[1]), but in the data you only provide one value init = c(1.0). I think BUGS gets confused here with the index (i.e. it expects two or more init values in the data, where you only provide one). The problem goes away when you replace init[1] with init in the model.

#BUGS model
model{
  solution[1:ngrid, 1] <-  ode(init,  tgrid[1:ngrid],    D(C[1], t),    origin,     tol)
  D(C[1], t) <- -lambda * C[1] 
  log.lambda ~ dnorm(0.0, tau.lambda); 
  lambda <- exp(log.lambda)
  for (i in 1:ngrid) {
    obs[i] ~ dnorm(solution[i, 1], tau) 
  }
  tau ~ dgamma(a, b)
}

#data
list(obs = c(0.7887, 1.241, 0.7051, 0.7388, 0.3903, 0.2632, 0.1174, 0.549, -0.1767, 0.02938, 0.154, 0.1465, 0.07878, -0.5569, 0.01293, 0.2905, -0.2665, -0.3881, 0.02517, 0.138, 0.4004, 0.2859, -0.1217, 0.3961, 0.3813, 0.1846, -0.3581, 0.3293, 0.04089, 0.01972, 0.3203, 0.5294, -0.1389, -0.3732, 0.1341, -0.02432, 0.2261, -0.3612, 0.3131, -0.258, 0.02948, -0.0208, 0.1066, 0.3796, -0.2645, 0.1035, 0.1001, -0.2415, 0.06739, -0.1554, -0.2388), 
tgrid = c(0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2,  2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6,  4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2,  7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8,  10),
origin = 0, tol = 0.001, ngrid = 51, init = 1, a = 0.001,  b = 0.001, tau.lambda = 0.01))

#inits
list(log.lambda = 1.0,tau = 0.01)

If you replace the model above with your previous model (Modele) everything should run from R smoothly.

General Point: I find it much easier to de-bug in OpenBUGS directly (using the script above and the GUI to compile, load data and inits) than via R2OpenBUGS.