[R-sig-eco] Bayesian Methods for Ecology: inits working with WinBUGS not working with R & R2WinBUGS
Guillaume Chapron
carnivorescience at gmail.com
Sun Dec 28 22:51:10 CET 2008
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
More information about the R-sig-ecology
mailing list