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
Bayesian Methods for Ecology: inits working with WinBUGS not working with R & R2WinBUGS
2 messages · Guillaume Chapron
Hello, The issue raised in my yesterday email is now solved. Philip Dixon emailed a suggestion to me first and I tested it and it works. Hereafter, you will find explanations, they are probably of interest to people switching between R and WinBUGS.
Remember that R and (Win)BUGS store matrices (and arrays) in different order. In BUGS, you specify values going across the rows. In R, you specify values going down the columns. If you have a reasonable matrix in R format, R2WinBUGS takes care of the ordering by transposing the matrix before passing it to WinBUGS. I suggest you print out your data and init matrices in R and check that they are in the order you intend. If they are in the wrong order, your NA's will NOT be where you intended, which explains the 'some variables uninitialized' message. -- So when pasting into R some WinBUGS code, it is important to remember this ordering issue. What I did is replacing this 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)) with this N = matrix(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), nrow = 4, ncol = 12, byrow = T) Cheers Guillaume