Hi everyone,

I recently did an analysis using WinBUGS to describe the vegetation
associated with a number of South American birds (Journal of Field
Ornithology 84:1-12).

Now I would like to do the same analysis but compare the vegetation where
birds forage with random sites in forest fragments where these birds are
(mostly) absent. I think it would be useful if I could show density plots
of the posterior estimates and potential overlap between bird habitat and
habitat in fragments.

I tried adapting my functional WinBUGS code to R2OpenBUGS without much
success. Here's the error that pops up in OpenBUGS:

model is syntactically correct
data loaded
model compiled
initial values loaded but chain contains uninitialized variables
error for node smtree[1] of type GraphPoisson.Node first argument must be
positive
model must be initialized before updating
model must be initialized before monitors used
model must be initialized before monitors used
model must be initialized before monitors used
model must be initialized before DIC can be monitored
model must be initialized before updating
model must be initialized before monitors used
DIC monitor not set

I'm using R 3.0.0 and OpenBUGS 3.2.2

The data and code are below. If anything can be clarified please let me
know.

Many thanks,

Jeff

structure(list(name = structure(c(1L, 1L, 1L, 1L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 8L,
8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L), .Label = c("CA", "con", "CT", "FA", "FC", "GV", "HFR", "HM",
"MC", "MF", "MT", "OFR", "TFR"), class = "factor"), site = c(1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L
), smtree = c(116L, 141L, 110L, 85L, 104L, 84L, 118L, 108L, 89L,
71L, 86L, 124L, 113L, 116L, 139L, 87L, 68L, 80L, 70L, 138L, 81L,
80L, 111L, 101L, 70L, 99L, 105L, 126L, 121L, 127L, 68L, 71L,
90L, 109L, 101L, 123L, 87L, 117L, 63L, 102L, 96L, 70L, 86L, 77L,
92L, 63L, 88L, 120L, 133L, 107L, 114L, 77L, 97L, 82L, 110L, 108L,
101L, 84L, 144L, 206L, 131L, 296L, 182L, 69L, 69L, 59L, 153L,
108L, 31L, 71L, 53L, 86L, 80L, 133L, 104L, 102L, 70L, 70L, 96L,
86L, 96L, 85L, 71L, 72L, 84L, 112L, 101L, 131L, 61L, 44L, 71L,
79L, 79L, 108L, 87L, 108L, 126L, 65L, 109L, 80L, 80L, 107L, 93L,
101L, 93L, 77L, 89L, 100L, 128L, 129L, 67L, 126L, 180L, 180L,
90L, 78L, 125L, 100L, 109L, 110L, 100L, 126L, 105L, 101L, 88L,
96L, 130L, 107L, 126L, 80L, 121L, 107L, 77L, 109L, 94L, 101L,
175L, 121L, 135L, 112L, 123L, 106L, 153L, 94L, 87L, 109L, 67L,
94L, 113L, 96L, 110L, 109L, 136L, 109L, 103L, 167L, 176L, 146L,
158L, 169L, 200L, 102L, 118L, 148L, 59L, 32L, 87L, 93L, 84L,
113L, 58L, 84L, 121L, 119L, 124L, 72L, 129L, 183L, 95L, 61L,
77L, 73L, 84L, 106L, 89L, 129L, 114L, 191L, 154L, 156L, 210L,
74L, 213L, 195L, 67L, 105L, 100L, 172L, 144L, 160L, 64L, 77L,
129L, 66L, 39L, 103L, 102L, 98L, 53L, 133L, 94L, 99L, 113L, 195L,
167L, 125L, 207L, 189L, 200L, 146L, 58L, 198L, 207L, 165L, 101L,
113L, 101L)), .Names = c("name", "site", "smtree"), class = "data.frame",
row.names = c(NA,
-227L))


# Write model
sink("anova.txt")
cat("

model{
for (i in 1:227)
{
re[i] ~ dnorm(0, tau) # extra-Poisson variation for data point i
mean[i] <- alpha[site[i]] + re[i] # log(mean density) = group mean+ extra
smtree[i] ~ dpois(mean[i]) # data drawn from Poisson
}
for (i in 1:13) # for each of the 12 groups
{
alpha[i] ~ dlnorm(0, 1.0E-6) # mean of group i (log scale)
median[i] <- exp(alpha[i])
}
sd ~ dunif(0, 10) # sd of the extra-Poisson variation
tau <- 1 / (sd*sd) # precision
}

",fill=TRUE)
sink()

inits <- function(){ list(alpha=c(0,0,0,0,0,0,0,0,0,0,0,0,0), sd=1)}

params <- c("median", "sd")

# MCMC settings
ni <- 10000
nb <- 2000
nt <- 2
nc <- 2

win.data <- list("site", "smtree")

# Start Gibbs sampling
out <- bugs(win.data, inits, params, "anova.txt", n.thin=nt, n.chains=1,
n.burnin=nb, n.iter=ni, debug = TRUE)

# Inspect estimates
print(out, dig = 3)


-- 
********************************************************
Jeffrey A. Stratford, PhD
Department of Biology and Health Sciences
Wilkes University, PA 18766 USA
570-332-2942
http://web.wilkes.edu/jeffrey.stratford/
********************************************************

	[[alternative HTML version deleted]]

