#24 Logic Error

[adigg]

I originally posted this on the discussion, but after much more testing think that it is a bug. I've since tried this for several other similar models, all of which run in BUGS, but give the same error in JAGS.

If anyone has any suggestions they'd be much appreciated.

Thanks.

---

JAGS gives the following error with a model that runs fine in OpenBUGS (via R2OpenBUGS):

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  LOGIC ERROR:
Invalid mask in LogicalNode::isClosed
Please send a bug report to martyn_plummer@users.sourceforge.net

(The email bounced.)

The model fits an AR1 correlation structure to a GAM, using an example from Zuur et al.'s 2012 book, "Zero-inflated models and generalized mixed models with R". When I try a Poisson response, it works OK, but JAGS produces the above error when I change to a binary response.

OS X 10.8.2
Jags 3.1.0
R 2.15.1
R2Jags 0.03-08

~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Jags example: GAM AR-1

From Zuur et al. 2012: "Zero Inflated Models and Generalized Linear Mixed Models with R", Chapter 9.8.2

Runs in OpenBUGS

require(R2jags)

Read data

SW <- read.table("http://www.highstat.com/Book4/Spermwhales.txt", header=T, sep="\t")
MyVar <- c("Year", "Nsea", "AnnualTTAnom", "WinterNAO")
SW1 <- SW[SW$Year >=1563 & SW$Year <= 2000, MyVar]
SW1<-SW1[1:dim(SW1)[1]-1,MyVar] # take off last row
SW1$TempAn <- as.numeric(scale(SW1$AnnualTTAnom))
SW1$NAO <- as.numeric(scale(SW1$WinterNAO))
N <- nrow(SW1)
num.knots <- 9
Prob <- seq(0, 1, length= (num.knots+2))
Prob1 <- Prob[-c(1, (num.knots+2))]
X <- cbind(rep(1,N), SW1$TempAn)
Knots <- quantile(unique(SW1$TempAn), probs=Prob1, na.rm=F)
GetZ <- function(x, Knots){
Z_K <- (abs(outer(x, Knots, "-")))^3
OMEGA_all <- (abs(outer(Knots, Knots, "-")))^3
svd.OMEGA_all <- svd(OMEGA_all)
sqrt.OMEGA_all <- t(svd.OMEGA_all$v %*% (t(svd.OMEGA_all$u) * sqrt(svd.OMEGA_all$d)))
Z <- t(solve(sqrt.OMEGA_all, t(Z_K)))
return(Z)
}
Z <- GetZ(SW1$TempAn, Knots)

win.data <- list(Y = SW1$Nsea, NumKnots=9, X=X, Z=Z, N=nrow(SW1))

Poisson GLMM with auto-regressive correlation:

Works

sink("zuur2012.9-8-2.pois.bugs.txt")
cat("
model{
# Priors for regression parameters
for (i in 1:2) {beta[i] ~ dnorm(0, 0.001)}
for (i in 1:NumKnots) {b[i] ~ dnorm(0, taub)}
taub ~ dgamma(0.01, 0.01)
taueps ~ dgamma(0.01, 0.01)
rho ~ dunif(-0.99, 0.99)

  # First row (i=1)
  # Calculate X*beta
  eta[1] <- beta[1] * X[1,1] + Smooth1[1]
  Smooth1[1] <- beta[2] * X[1,2] + b[1]*Z[1,1] + b[2]*Z[1,2] + b[3]*Z[1,3] + b[4]*Z[1,4] + b[5]*Z[1,5] + b[6]*Z[1,6] + b[7]*Z[1,7] + b[8]*Z[1,8] + b[9]*Z[1,9]
  # Calculate X*beta + eps 
  eta1[1] ~ dnorm(eta[1], taueps)
  log(mu[1]) <- eta1[1]
  Y[1] ~ dpois(mu[1])
  Pres[1] <- (Y[1] - mu[1])/pow(mu[1], 0.5)  # Pearson residuals

  for(i in 2:N){
  # Binomial part
  eta[i] <- beta[1] * X[i,1] + Smooth1[i]
  Smooth1[i] <- beta[2] * X[i,2] + b[1]*Z[i,1] + b[2]*Z[i,2] + b[3]*Z[i,3] + b[4]*Z[i,4] + b[5]*Z[i,5] + b[6]*Z[i,6] + b[7]*Z[i,7] + b[8]*Z[i,8] + b[9]*Z[i,9]
    # Calculate eta + rho*(log(u[i=1]) - eta[i-1]) + eps
    temp[i] <- eta[i] + rho * (log(mu[i-1]) - eta[i-1])
    eta1[i] ~ dnorm(temp[i], taueps)
    log(mu[i]) <- max(-10, min(10, eta1[i]))
    Y[i] ~ dpois(mu[i])
    Pres[i] <- (Y[i] - mu[i])/pow(mu[i], 0.5)  # Poisson residuals

  }
}
", fill=TRUE)

sink()

inits <- function(){
list( beta = rnorm(2, 0, 0.001), b = rnorm(9, 0, 0.001), taub=runif(1, 0.1, 0.5), taueps=runif(1, 0.1, 5), rho=runif(1, -0.99, 0.99))
}

params <- c("beta", "b", "Smooth1", "rho", "Pres")

nb<-500
nit<-3000
nt<-10
mjags.pois <- jags(data=win.data, inits=inits, parameters.to.save=params, model.file="zuur2012.9-8-2.pois.bugs.txt", n.chains=3, n.burnin=nb, n.iter=nit, n.thin=nt)

Same, but with binary response variable

Crashes with LOGIC ERROR

win.data2<-win.data
win.data2$Y[win.data2$Y > 0] <-1 # Make response binary
sink("zuur2012.9-8-2.binar.bugs.txt")
cat("
model{
# Priors for regression parameters
for (i in 1:2) {beta[i] ~ dnorm(0, 0.001)}
for (i in 1:NumKnots) {b[i] ~ dnorm(0, taub)}
taub ~ dgamma(0.01, 0.01)
taueps ~ dgamma(0.01, 0.01)
rho ~ dunif(-0.99, 0.99)

  # First row (i=1)
  # Calculate X*beta
  eta[1] <- beta[1] * X[1,1] + Smooth1[1]
  Smooth1[1] <- beta[2] * X[1,2] + b[1]*Z[1,1] + b[2]*Z[1,2] + b[3]*Z[1,3] + b[4]*Z[1,4] + b[5]*Z[1,5] + b[6]*Z[1,6] + b[7]*Z[1,7] + b[8]*Z[1,8] + b[9]*Z[1,9]
  # Calculate X*beta + eps 
  eta1[1] ~ dnorm(eta[1], taueps)
  logit(mu[1]) <- eta1[1]
  Y[1] ~ dbin(mu[1],1)

  for(i in 2:N){
  eta[i] <- beta[1] * X[i,1] + Smooth1[i]
  Smooth1[i] <- beta[2] * X[i,2] + b[1]*Z[i,1] + b[2]*Z[i,2] + b[3]*Z[i,3] + b[4]*Z[i,4] + b[5]*Z[i,5] + b[6]*Z[i,6] + b[7]*Z[i,7] + b[8]*Z[i,8] + b[9]*Z[i,9]
    # Calculate eta + rho*(logit(u[i=1]) - eta[i-1]) + eps
    temp[i] <- eta[i] + rho * (logit(mu[i-1]) - eta[i-1])
    eta1[i] ~ dnorm(temp[i], taueps)
    logit(mu[i]) <- max(-10, min(10, eta1[i]))
    Y[i] ~ dbin(mu[i],1)

  }
}
", fill=TRUE)

sink()

inits <- function(){
list( beta = rnorm(2, 0, 0.001), b = rnorm(9, 0, 0.001), taub=runif(1, 0.1, 0.5), taueps=runif(1, 0.1, 5), rho=runif(1, -0.99, 0.99))
}
params <- c("beta", "b", "Smooth1", "rho")

nb<-500
nit<-3000
nt<-10

mjags.binar <- jags(data=win.data2, inits=inits, parameters.to.save=params, model.file="zuur2012.9-8-2.binar.bugs.txt", n.chains=3, n.burnin=nb, n.iter=nit, n.thin=nt, DIC=T)

~~~~~~~~~~~~~~~~~~~~~~~~