
Hi,

I'm not sure if this is the right forum to ask questions about MCMC; apologies if not.

I'm using R2jags and want to calculate Pearson residuals individually for each MCMC realisation. However, when I use the appropriate WinBUGS code to do this in JAGS using R2jags, I get an error saying that the node for the Pearson residual variable is not found:

> Error in jags.samples(model, variable.names, n.iter, thin, type = "trace",  : 
>  Failed to set trace monitor for PRes
>  Node not found

However, I have defined PRes in the model file. My model has two levels of nesting, which I think is the source of the problem: with a single nest the residuals are calculated OK (e.g. the example R2jags example at http://people.ucsc.edu/~abrsvn/bayes_winbugs_jags_5.r works).

I'd be very grateful for any suggestions on how to fix this. I'm new to MCMC/jags, so please excuse me if the answer's obvious. I'm using OS X and can't get WinBUGs working, so would prefer to stick with JAGS.

Below is a reproducible example, copied from Chapter 2 Zuur et al.'s 2012 book, "Zero-inflated models and generalized linear mixed models with R". (Apologies for the rather long code, and to Alain Zuur et al. for using their code.) The example works without 'P.Res' in the parameters; adding it causes the jags() call to fail. The offending line in the model is:

      PRes[j,i] <- (NCalls[j,i] - mu[j,i]) / sqrt(mu[j,i])  


Thanks very much,

Andrew





# Zuur et al. 2012 Ch2: Owl data

# Owl data
Owls <- read.table("http://www.highstat.com/Book4/Owls.txt", header=T)
Owls$NCalls <- Owls$SiblingNegotiation

# Remove nests with <=10 obs:
Owls2 <- Owls[Owls$Nest != "Chevroux" & Owls$Nest != "Forel" & Owls$Nest != "GDLV" & Owls$Nest != "SEvaz", ]
Owls2$Nest <- factor(Owls2$Nest)
Owls2$LogBroodSize <- log(Owls2$BroodSize)

#  2.6.1
# Convert variable vectors to matrices, with each column a nest, and rows the observations for each

Vec2Mat <- function(FNest, FY){
  AllNests <- levels(FNest)
  NNest <- length(unique(FNest))
  MaxObs <- max(table(FNest))
  Y.ij <- matrix(NA, nrow=MaxObs, ncol=NNest)
  for (i in 1:NNest){
    Selection <- FNest == AllNests[i]
    ni <- sum(Selection)
    Y.ij[1:ni, i] <- FY[Selection]
  }
  Y.ij  
}

NCalls.ij <- Vec2Mat(Owls2$Nest, Owls2$NCalls)
LBroodSize.ij <- Vec2Mat(Owls2$Nest, Owls2$LogBroodSize)
SexParent.ij <- Vec2Mat(Owls2$Nest, Owls2$SexParent)
FoodTreatment.ij <- Vec2Mat(Owls2$Nest, Owls2$FoodTreatment)

# Centre cts variables:
Owls2$Cen.AT <- Owls2$ArrivalTime - mean(Owls2$ArrivalTime)
ArrivalTime.ij <- Vec2Mat(Owls2$Nest, Owls2$Cen.AT)

# Categorical variables as 0/1:
iSexParent.ij <- SexParent.ij - 1
iFoodTreatment.ij <- FoodTreatment.ij - 1

# Number obs per nest:
NObservationsInNest <- as.numeric(tapply(Owls2$NCalls, FUN=length, INDEX=Owls2$Nest))


#####  2.6.2 WinBUGS GLM
owldata <- list(NCalls = NCalls.ij, ArrivalTime = ArrivalTime.ij, LBroodSize = LBroodSize.ij, NNest = length(levels(Owls2$Nest)), NObservationsInNest = NObservationsInNest, iSexParent = iSexParent.ij, iFoodTreatment=iFoodTreatment.ij)

# Model:
sink("modelglm.txt")
cat("
model{
  
  #Priors
  for (i in 1:5) {beta[i] ~ dnorm(0, 0.001)}
  alpha ~ dnorm(0, 0.001)
  
  # Likelihood
  for (i in 1:NNest){
    for (j in 1:NObservationsInNest[i]){
      NCalls[j,i] ~ dpois(mu[j,i])
      log(mu[j,i]) <- max(-20, min(20, eta[j,i]))
      eta[j,i] <- alpha + beta[1]*iSexParent[j,i] + beta[2]*iFoodTreatment[j,i] + beta[3]*ArrivalTime[j,i] + beta[4]*iFoodTreatment[j,i]*iSexParent[j,i] + beta[5]*iSexParent[j,i]*ArrivalTime[j,i] + 1*LBroodSize[j,i]

      YNew[j,i] ~ dpois(mu[j,i])  
      PRes[j,i] <- (NCalls[j,i] - mu[j,i]) / sqrt(mu[j,i])  
      PResNew[j,i] <- (YNew[j,i] - mu[j,i]) / sqrt(mu[j,i]) # Pearson resids for replicate dataset
      D[j,i] <- pow(PRes[j,i], 2)
      DNew[j,i] <- pow(PResNew[j,i], 2)
    }
  Fiti[i] <- sum(D[1:NObservationsInNest[i], i])       # Add up discrepancy measures for each obs
  FitiNew[i] <- sum(DNew[1:NObservationsInNest[i], i])
  }
  Fit <- sum(Fiti[1:NNest])          # Add up discrepancy measures for each nest
  FitNew <- sum(FitiNew[1:NNest])
}
", fill=T)
sink()

# Initialise chains:
jags.inits <- function() { list(alpha=rnorm(1), beta = rnorm(5))}
jags.params <- c("alpha", "beta", "Fit", "FitNew")
jags.params <- c("alpha", "beta", "PRes","Fit", "FitNew")

# Run JAGS model:
require(R2jags)
jags.glm <- jags(data=owldata, inits=jags.inits, parameters.to.save=jags.params, model="modelglm.txt", n.thin=10, n.chains=3, n.burnin=100, n.iter=1000)

	[[alternative HTML version deleted]]

