1

I am using R2 Winbugs to run an open-population binomial mixture model following Kery et al 2009 (paper here) with real data on surveys.

As I first try I am using only one covariate for the abundance (X1) that I made up myself. The model compiles and runs but I get "undefined real results".

Could anyone help me at finding how to fix this if possible?

My count data can be downloaded here

The code is as follows

#data formatting
bdat <- read.table("PyroniaC_model_test2.txt", header = TRUE)
names(bdat)
#convert data to array
y <- array(NA, dim = c(620, 2, 352))  # 620 sites, 2 reps, 352 days
str(y)
 for(k in 1:352){
 sel.rows <- bdat$Bweekcoded == k
 y[ , , k] <- as.matrix((bdat)[sel.rows, 3:4],stringsAsFactors = F)
 }
#covs
X1<-rnorm(620, 0,100)

# Specify model in BUGS language
sink("Nmix3C.txt")
cat("
model {

# Priors
alpha0 ~ dunif(0, 0.2) #intercept log(N)
beta1 ~ dunif(0, 2) #beta cov1
r.rate ~ dunif(-1, 1) #growth rate
for (k in 1:352){   #prior for detection
 p[k] ~ dunif(0, 1)
}
# Likelihood
# Ecological model for true abundance
for (k in 1:352){                          # Loop over weeks (352)
for (i in 1:R){                       # Loop over R sites (620)
log(lambda[i,k])<-alpha0+beta1*X1[i]+r.rate*(k-1)
N[i,k] ~ dpois(lambda[i,k])          # Abundance

# Observation model for replicated counts
for (j in 1:T){                    # Loop over temporal reps (2)
y[i,j,k] ~ dbin(p[k], N[i,k])   # Detection

# Assess model fit using Chi-squared discrepancy
# Compute fit statistic E for observed data
eval[i,j,k] <- p[k] * N[i,k]   # Expected values
E[i,j,k] <- pow((y[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k] + 0.5)
# Generate replicate data and compute fit stats for them
y.new[i,j,k] ~ dbin(p[k], N[i,k])
E.new[i,j,k] <- pow((y.new[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k] + 0.5)

} #j 
} #i 
} #k 

# Derived and other quantities
for (k in 1:352){
totalN[k] <- sum(N[,k])# Total pop. size across all sites
#mean.abundance[k] <- exp(alpha.lam[k])
}
fit <- sum(E[,,])
fit.new <- sum(E.new[,,])
}

",fill = TRUE)
sink()

# Bundle data
R = nrow(y)
T = ncol(y)
win.data <- list(y = y, R = R, T = T, X1=X1)

 # Initial values
Nst <- apply(y, c(1, 3), max) +1
Nst[is.na(Nst)] <- 1
inits <- function(){list(N = Nst , alpha0 = runif(1, 0, 0.2),r.rate=runif(1,-1,1),beta1 = runif(1,     0,2), p=runif(1,0,1))}

# Parameters monitored
params <- c("totalN", "alpha0", "beta1","p", "fit", "fit.new", "r.rate")

# MCMC settings
ni <- 10000
nt <- 8
nb <- 2000
nc <- 3

library(R2WinBUGS)
# Call WinBUGS from R (BRT 1 min)
out3C <- bugs(win.data, inits, params, "Nmix3C.txt", n.chains = nc, n.thin = nt, n.iter = ni,  n.burnin = nb, debug = TRUE)

Many thanks in advance to all of you!

PS. I also tried with a real covariate X1 but I get the same error.

YMC
  • 63
  • 6
  • Have you tried some of the some of the suggestions here... http://stackoverflow.com/questions/21969600/bugs-error-messages – guyabel Jan 14 '15 at 08:33
  • Yes, but it did not work. I have started from zero with the code helped by a colleague of mine and it seems the new one is working. Apart of the fact that there were some problem with the counts (needed to be transformed to integers)I still don't know why it does not work. The strange thing is that even copying the code from the paper I get the same error in WINBUGS but it works on JAGS. Thanks for answering, btw!! – YMC Jan 15 '15 at 09:04

0 Answers0