[R-sig-eco] Bayesian Methods for Ecology: inits working with WinBUGS not working with R & R2WinBUGS

[Guillaume Chapron]

Dear all,

I'm trying to have a WinBUGS example from the book "Bayesian Methods  
for Ecology" by Michael A. McCarthy running with R through R2WinBUGS.  
I cannot have this working, as there seems to be an issue with inits.

The model is a population dynamic one from the book Chapter 9, in  
particular Box 9.1 p.219. The WinBUGS code is as follows:

##################################
# BUGS model to be called from WinBUGS
##################################

model
{
   for (j in 1:4)  # for each population
   {
     for(i in 1:T[j])  # for each year
     {
# env stoch in population growth rate, drawn from normal
       ev[j, i] ~ dnorm(0, tau)
# per capita growth rate is dens dep - Ricker
       Er[j, i] <- exp(alpha*(1 - N[j, i]/K[j]) + ev[j, i])
# lambda equal to number this year times per capita rate
       lambda[j, i] <- N[j, i] * Er[j, i]
# number next year drawn from Poisson with mean lambda - demo. stoch.
       N[j, i+1] ~ dpois(lambda[j, i])
     }
  }

#PRIORS
# alpha is maximum exponential growth rate	
   alpha ~ dunif(0, 1.0986)

# K's are carrying capacities of the 4 sites	
   for (i in 1:4)
   {
     K[i] ~ dunif(1, 50)
   }

# st. dev. in growth rate due to env. stoch.
   sd ~ dunif(0, 0.5)
   tau <- 1/sd/sd
}


list(alpha = 0.5, K=c(30,30, 30, 30), sd=0.1, ev = structure(.Data = c(
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1 
,            NA,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,  NA, 
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,  NA),
.Dim = c(4,11)))


list(N = structure(.Data = c(
32, 28, 29, 39, 20, 24, 22, 35, 20, 36, 34, 40,
25, 25, 24, 24, 28, 15, 30, 36, 27, 31, 20, NA,
7, 11, 8, 14, 11, 5, 12, 6, 10, 12, 19, NA,
12, 16, 8, 8, 16, 11, 10, 6, 4, 10, 12, NA),
.Dim = c(4, 12)), T = c(11, 10, 10, 10))

##################################

I have explicitly written the inits in this case, I have just rounded  
the ones given by WinBUGS. All this works perfectly well.

Running this with R would be done by this code (if I'm correct):

##################################
# R script to call R2WinBUGS and plot density
##################################

# Define the dataset

mydata <- list(N = structure(.Data = c(32, 28, 29, 39, 20, 24, 22, 35,  
20, 36, 34, 40, 25, 25, 24, 24, 28, 15, 30, 36, 27, 31, 20, NA,7, 11,  
8, 14, 11, 5, 12, 6, 10, 12, 19, NA, 12, 16, 8, 8, 16, 11, 10, 6, 4,  
10, 12, NA),.Dim = c(4, 12)), T = c(11, 10, 10, 10))

# Define inits

inits <- list(list(alpha = 0.5, K=c(30,30, 30, 30), sd=0.1, ev =  
structure(.Data =  
c 
(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1 
,NA,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,NA, 
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,NA), .Dim = c(4,11))))

# specify the parameters to be monitored
parameters <- c("alpha", "K", "sd")

# load R2WinBUGS package
library(R2WinBUGS)

# run the MCMC analysis
res.sim <- bugs(data=mydata, inits, parameters,"pva.bug",  
n.chains=length(inits), n.burnin=5000, n.iter=10000
#, bugs.directory="/Applications/WinBUGS14", working.directory=".",  
#WINE="/Applications/Darwine/Wine.bundle/Contents/bin/wine",  
#WINEPATH="/Applications/Darwine/Wine.bundle/Contents/bin/winepath"
)

attach.bugs(res.sim)
plot(density(alpha))

##################################

and then the BUG model named pva.bug would be this (same as before  
without the inits and data in fact)

##################################
# BUGS model to be called from R
##################################

model
{
   for (j in 1:4)  # for each population
   {
     for(i in 1:T[j])  # for each year
     {
# env stoch in population growth rate, drawn from normal
       ev[j, i] ~ dnorm(0, tau)
# per capita growth rate is dens dep - Ricker
       Er[j, i] <- exp(alpha*(1 - N[j, i]/K[j]) + ev[j, i])
# lambda equal to number this year times per capita rate
       lambda[j, i] <- N[j, i] * Er[j, i]
# number next year drawn from Poisson with mean lambda - demo. stoch.
       N[j, i+1] ~ dpois(lambda[j, i])
     }
  }

#PRIORS
# alpha is maximum exponential growth rate
   alpha ~ dunif(0, 1.0986)

# K's are carrying capacities of the 4 sites
   for (i in 1:4)
   {
     K[i] ~ dunif(1, 50)
   }

# st. dev. in growth rate due to env. stoch.
   sd ~ dunif(0, 0.5)
   tau <- 1/(sd*sd)

}

##################################

BUT, this does not work, the log file says "this chain contains  
uninitialized variables" and later "cannot bracket slice for node  
N[4,6]". I cannot understand why something perfectly working directly  
in WinBUGS (and giving the correct results) will not work at all when  
WinBUGS is called by R2WinBUGS. What am I doing wrong?

Thanks for your help !
Cheers

Guillaume