Error OpenBUGS undefined

Question

I am working on a binomial mix model using the OpenBUGS and R package R2OpenBUGS. I have successfully created simpler models, but as soon as I add another level for imperfect detection, I consistently get an error variable X is not defined in model or in data set. I tried several different things, including changing the structure of my data and entering my data directly into OpenBUGS. I post this in the hope that someone else has experience with this error and may know why OpenBUGS does not recognize the variable X, even if it is clearly defined, as far as I can tell.

I also got an error expected the collection operator c error pos 8- this is not an error I received earlier, but I am also confused.

Both the model and the data modeling function are taken from Kery Introduction to WinBUGS for Ecologists (2010). I note that the data set is here instead of my own data, which is similar.

I turn on the function of creating a dataset as well as a model. Sorry for the length.

# Simulate data: 200 sites, 3 sampling rounds, 3 factors of the level 'trt', 
# and continuous covariate 'X'

data.fn <- function(nsite = 180, nrep = 3, xmin = -1, xmax = 1, alpha.vec = c(0.01,0.2,0.4,1.1,0.01,0.2), beta0 = 1, beta1 = -1, ntrt = 3){
  y <- array(dim = c(nsite, nrep))  # Array for counts
  X <- sort(runif(n = nsite, min = xmin, max = xmax))   # covariate values, sorted
  # Relationship expected abundance - covariate
  x2 <- rep(1:ntrt, rep(60, ntrt)) # Indicator for population
  trt <- factor(x2, labels = c("CT", "CM", "CC"))
  Xmat <- model.matrix(~ trt*X)
  lin.pred <- Xmat[,] %*% alpha.vec # Value of lin.predictor
  lam <- exp(lin.pred)
  # Add Poisson noise: draw N from Poisson(lambda)
  N <- rpois(n = nsite, lambda = lam)
  table(N)                # Distribution of abundances across sites
  sum(N > 0) / nsite          # Empirical occupancy
  totalN <- sum(N)  ;  totalN
  # Observation process
  # Relationship detection prob - covariate
  p <- plogis(beta0 + beta1 * X)
  # Make a 'census' (i.e., go out and count things)
  for (i in 1:nrep){
    y[,i] <- rbinom(n = nsite, size = N, prob = p)
  }
  # Return stuff
  return(list(nsite = nsite, nrep = nrep, ntrt = ntrt, X = X, alpha.vec = alpha.vec, beta0 = beta0, beta1 = beta1, lam = lam, N = N, totalN = totalN, p = p, y = y, trt = trt))
}

data <- data.fn()

And here is the model:

sink("nmix1.txt")
cat("
    model {

    # Priors
    for (i in 1:3){     # 3 treatment levels (factor)   
    alpha0[i] ~ dnorm(0, 0.01)       
    alpha1[i] ~ dnorm(0, 0.01)       
    }
    beta0 ~ dnorm(0, 0.01)       
    beta1 ~ dnorm(0, 0.01)

    # Likelihood
    for (i in 1:180) {      # 180 sites
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha0[trt[i]] + alpha1[trt[i]]*X[i]

    for (j in 1:3){     # each site sampled 3 times
    y[i,j] ~ dbin(p[i,j], C[i])
    lp[i,j] <- beta0 + beta1*X[i]
    p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
    }
    }

    # Derived quantities

    }
    ",fill=TRUE)
sink()

# Bundle data
trt <- data$trt
y <- data$y
X <- data$X
ntrt <- 3

# Standardise covariates
s.X <- (X - mean(X))/sd(X)

win.data <- list(C = y, trt = as.numeric(trt), X = s.X)

# Inits function
inits <- function(){ list(alpha0 = rnorm(ntrt, 0, 2), 
                          alpha1 = rnorm(ntrt, 0, 2),
                beta0 = rnorm(1,0,2), beta1 = rnorm(1,0,2))}

# Parameters to estimate
parameters <- c("alpha0", "alpha1", "beta0", "beta1")

# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3

# Start Markov chains
out <- bugs(data = win.data, inits, parameters, "nmix1.txt", n.thin=nt, 
            n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

Answer 1

Note. . This answer went through a major revision after I noticed another problem with the code.

---

, y N C . y () C Bugs, . , , C "" ( ), .. N . y () , Bugs.

, , :

sink("nmix1.txt")
cat("
    model {

    # Priors
    for (i in 1:3){     # 3 treatment levels (factor)   
    alpha0[i] ~ dnorm(0, 0.01)       
    alpha1[i] ~ dnorm(0, 0.01)       
    }
    beta0 ~ dnorm(0, 0.01)       
    beta1 ~ dnorm(0, 0.01)

    # Likelihood
    for (i in 1:180) {      # 180 sites
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha0[trt[i]] + alpha1[trt[i]]*X[i]

    for (j in 1:3){     # each site sampled 3 times
        y[i,j] ~ dbin(p[i,j], C[i])
        lp[i,j] <- beta0 + beta1*X[i]
        p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
    }
    }

    # Derived quantities

    }
    ",fill=TRUE)
sink()

# Bundle data
trt <- data$trt
y <- data$y
X <- data$X
N<- data$N
ntrt <- 3

# Standardise covariates
s.X <- (X - mean(X))/sd(X)

win.data <- list(y = y, trt = as.numeric(trt), X = s.X, C= N)

# Inits function
inits <- function(){ list(alpha0 = rnorm(ntrt, 0, 2), 
                          alpha1 = rnorm(ntrt, 0, 2),
                beta0 = rnorm(1,0,2), beta1 = rnorm(1,0,2))}

# Parameters to estimate
parameters <- c("alpha0", "alpha1", "beta0", "beta1")

# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3

# Start Markov chains
out <- bugs(data = win.data, inits, parameters, "nmix1.txt", n.thin=nt, 
            n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

, - -1. 1 (~ = -0.4), Bugs, (.. ). , , , .

Answer 2

, OpenBUGS. ,

Ndata <- list(yrs=N$yrs, site=N$site), ... )

"" "". OpenBUGS

, ,

Ndata <- list(yrs=N$yrs, site=as.numeric(N$site)), ... )