Jolly-Seber runs with WinBUGS not with Jags

Question

I can run a simple Jolly-Seber model with WinBUGS but not with Jags. I can run a linear regression with Jags, which suggests to me that R is able to locate and execute Jags. Therefore, I suspect the problem may be that Jags cannot interpret one (or more) lines in the model code. Please examine the code below and suggest how it can be modified to run in Jags. Initially I suspected perhaps the prod function is not available in Jags. However, a search of the Jags manual shows that Jags does include a prod function.

This is the simplest working example I can think of, but will simplify further if possible.

An example data set is provided at the bottom. The model code is modified from Kery and Schaub (2012).

# BUGS code

sink("C:/Users/mmiller/Documents/simple R programs/my.model.txt")
cat("
model {
for (i in 1:M) {
   for (t in 1:(n.occasions-1)) {
      phi[i,t] <- mean.phi
   }
   for (t in 1:n.occasions) {
      p[i,t] <- mean.p
   }
}
mean.phi ~ dunif(0, 1)
mean.p   ~ dunif(0, 1)
for (t in 1:n.occasions) {
   gamma[t] ~ dunif(0, 1)
}
for (i in 1:M) {
   z[i,1] ~ dbern(gamma[1])
   mu1[i] <- z[i,1] * p[i,1]
   y[i,1] ~ dbern(mu1[i])

   for (t in 2:n.occasions) {
      q[i,t-1] <- 1-z[i,t-1]
      mu2[i,t] <- phi[i,t-1] * z[i,t-1] + gamma[t] * prod(q[i,1:(t-1)])
      z[i,t] ~ dbern(mu2[i,t])
      mu3[i,t] <- z[i,t] * p[i,t]
      y[i,t] ~ dbern(mu3[i,t])
   }
}
}
",fill=TRUE)
sink()

# run R2WinBUGS

setwd('C:/Users/mmiller/Documents/simple R programs')

library(R2WinBUGS)
data <- list(y = my.data, n.occasions = dim(my.data)[2], M = dim(my.data)[1])
inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = my.data)}  
parameters <- c("mean.p", "mean.phi")
bugs.out <- bugs(data, inits, parameters, 
                 "C:/Users/mmiller/Documents/simple R programs/my.model.txt", 
                 n.thin=1, n.chains=2, n.burnin=500, n.iter=1000, debug=FALSE, working.dir=getwd())

print(bugs.out, digits=2)

# run R2jags

library('R2jags')
data <- list(y = my.data, n.occasions = dim(my.data)[2], M = dim(my.data)[1])
inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = my.data)}  
parameters <- c("mean.p", "mean.phi")

jags.out2 <- jags(data = data, inits = inits, parameters, 
                  model.file = "C:/Users/mmiller/Documents/simple R programs/my.model.txt",
                  n.thin=1, n.chains=2, n.burnin=500, n.iter=1000, working.dir=getwd())
print(jags.out2, digits=2)

Here is the example data set:

my.data <- read.table(text = '

    1    1    1    0    0    0    0
    0    1    1    1    1    0    0
    0    1    0    0    0    0    0
    0    0    1    1    0    1    0
    0    1    1    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    1    0    0    0    0
    1    0    0    0    0    0    0
    1    0    1    1    1    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    1    0    0    0    0    0
    0    0    0    0    1    0    0
    1    0    0    0    0    0    0
    1    1    0    0    0    0    0
    1    0    1    0    0    0    0
    0    1    0    0    0    0    0
    0    1    1    0    0    0    0
    1    0    0    0    0    0    0
    0    1    1    0    0    0    0
    0    0    1    0    0    0    0
    0    0    0    1    0    0    0
    0    1    0    0    0    0    0
    0    0    0    0    1    1    0
    0    1    1    0    0    0    0
    0    1    1    0    0    0    0
    0    0    1    1    0    0    0
    0    0    0    1    0    0    0
    0    0    1    1    0    0    0
    0    0    1    1    1    0    0
    0    0    1    1    0    0    0
    0    0    1    0    0    0    0
    0    0    0    1    0    1    1
    0    0    0    0    1    1    0
    0    0    0    1    0    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    1    1
    0    0    0    1    0    0    0
    0    0    0    1    0    1    1
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    0    0    1    1
    0    0    0    0    0    0    1
    0    0    0    0    1    1    0
    0    0    0    0    0    1    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    0
    0    0    0    0    0    1    0
    0    0    0    0    0    1    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    0
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1', header = FALSE)

nz <- 300
my.data <- rbind(my.data, matrix(0, ncol = ncol(my.data), nrow = nz))
dim(my.data)
head(my.data)
my.data <- as.matrix(my.data)

Here is a linear regression that does run in Jags:

# Linear regression in JAGS using R2jags

library('R2jags')

x <- rnorm(10)
mu <- -3.2 + 1.5 * x
y <- rnorm(10, mu, sd = 4)

cat ("
model {
  for (i in 1:10) {
     y[i] ~ dnorm(mu[i], tau)
     mu[i] <- beta0 + beta1*x[i]
   }

  beta0 ~ dnorm(0, .01)
  beta1 ~ dnorm(0, .01)
  sigma ~ dunif(0,100)
  tau <- 1 / (sigma * sigma)

}

", file = 'C:/Users/mmiller/Documents/simple R programs/normal.txt')

data <- list(y=y, x=x)

inits <- function()

list(beta1 = rnorm(1), beta0 = rnorm(1), sigma = runif(1,0,2))

parameters <- c("beta0", "beta1", "sigma", "tau")

out <- jags(data = data, inits = inits, parameters, 
            model.file = 'C:/Users/mmiller/Documents/simple R programs/normal.txt',
            n.thin=1, n.chains=2, n.burnin=2000, n.iter=6000, working.dir=getwd())

print(out, digits=2)

Jolly-Seber runs with WinBUGS not with Jags

1 Answers