JAGS code for estimating group means with Beta distribution

Question

I'd like to estimate the means and sd's of percent canopy cover for 13 sites (9 are birds and 4 are potential habitats) using JAGS. I'm using a beta distribution to account for the fact that the data are bound by 0 and 1.

I have code for the model statement that works perfectly for other distributions (Poisson and log-normal) and I was attempting to adapt that code but I failed miserably.

Below are the R code, the model statement, and the data. I'm using R 3.1.1 in Windows Vista. If you could look at the model statement and straighten me out I would be very thankful.

Thanks,

Jeff

######## MODEL ##############
model{
  for (i in 1:227) {
    log(mean[i]) <- a[site[i]] 
    cover20p[i] ~ dbeta(1, 0.5)   
  }
  for (i in 1:13){
    a[i] ~ dnorm(0, tau) 
    median[i] <- exp(a[i])
  }
  sd ~ dunif(0, 10) 
  tau <- 1 / (sd*sd) # precision
} 

#########  R  code ########## 
frag <- read.csv("f:\\brazil\\TIandFRAG.csv", header=T)
library(R2jags)
library(rjags)
setwd("f://brazil")
site <- frag$site
cover20p <- frag$cover20p/100
N <- length(frag$site)

jags.data <- list("site", "cover20p")
jags.params <- c("median", "test100MF","test100MT","test100fc","test100fa", 
"test100gv","test100hm","test100mc", "test100ca","test100ct", "test10MF",
"test10MT", "test10fc","test10fa", "test10gv", "test10hm", "test10mc", "test10ca", "test10ct", 
"test1MF", "test1MT", "test1fc",  "test1fa",  "test1gv", "test1hm", 
"test1mc", "test1ca", "test1ct", "t1est1_con","t2est10_con","t3est100_con",
"t4est1_100","t5est1_10","t6est10_100")
#inits1 <- list(a=0, sd=0)
#inits2 <- list(a=100, sd=50)
#jags.inits <- list(inits1, inits2)

jags.inits <- function() {
  list(a=c(0,0,0,0,0,0,0,0,0,0,0,0,0), sd=1)}

jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
n.iter=1000000, n.burnin=20000, model.file="fragmodelbeta.txt")

my.coda <- as.mcmc(jagsfit)
summary(my.coda, quantiles=c(0.05, 0.25,0.5,0.75, 0.95))
print(jagsfit, digits=3)

##### DATA ###################    
structure(list(site = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 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, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 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, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 
10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L
), canopy = c(0.95, 0.8, 0.85, 0.9, 0.35, 0.999, 0.999, 0.999, 
0.95, 0.55, 0.9, 0.85, 0.7, 0.65, 0.05, 0.6, 0.999, 0.999, 0.85, 
0.9, 1e-04, 0.45, 0.999, 0.7, 0.95, 0.5, 0.95, 0.6, 0.65, 0.7, 
0.4, 0.85, 0.6, 0.95, 0.75, 0.9, 0.85, 0.75, 0.7, 0.85, 0.3, 
0.7, 0.8, 0.7, 0.75, 0.8, 0.75, 0.95, 0.9, 0.05, 0.85, 0.6, 0.65, 
0.5, 0.85, 0.95, 0.85, 0.25, 0.75, 0.999, 0.65, 0.95, 0.8, 0.9, 
0.6, 0.8, 0.999, 0.2, 0.8, 0.4, 0.999, 0.95, 0.4, 0.999, 0.999, 
0.95, 0.45, 0.2, 0.7, 0.95, 0.7, 0.8, 0.5, 0.85, 0.55, 1e-04, 
0.25, 0.45, 0.999, 0.95, 0.999, 0.9, 0.6, 0.35, 0.95, 0.3, 0.999, 
0.999, 0.5, 0.4, 0.9, 0.999, 0.7, 0.999, 0.9, 0.999, 0.4, 0.55, 
0.8, 0.7, 0.999, 1e-04, 0.8, 1e-04, 0.7, 0.5, 0.8, 0.75, 1e-04, 
0.45, 0.1, 1e-04, 0.4, 0.55, 0.4, 0.999, 0.9, 0.9, 0.15, 0.55, 
0.35, 0.9, 0.65, 0.25, 0.999, 0.85, 0.999, 0.95, 0.7, 0.5, 0.7, 
0.2, 0.95, 0.999, 0.999, 0.25, 0.85, 0.5, 0.8, 0.75, 0.85, 0.7, 
0.95, 0.05, 0.65, 0.65, 0.999, 0.999, 0.999, 0.65, 0.4, 0.6, 
0.9, 0.85, 0.75, 0.5, 0.65, 0.999, 0.65, 0.55, 0.75, 0.4, 0.9, 
0.35, 0.999, 0.999, 0.4, 0.5, 0.8, 0.95, 0.95, 0.55, 0.7, 0.85, 
0.8, 0.8, 0.65, 0.999, 0.6, 0.5, 0.999, 0.8, 0.999, 0.45, 0.999, 
0.999, 0.8, 0.85, 0.999, 0.999, 0.999, 0.999, 0.5, 0.6, 0.15, 
0.75, 0.6, 0.1, 0.05, 1e-04, 0.999, 0.6, 0.1, 0.35, 0.9, 0.9, 
0.95, 0.95, 0.9, 0.55, 0.65, 0.9, 0.4, 0.999, 0.65, 0.5, 0.8)), .Names = c("site", 
"canopy"), class = "data.frame", row.names = c(NA, -227L))

Answer 1

In your model, you have cover20p as one of the variables, but have the data for canopy in the frag data.frame. I suspect you want canopy[i] ~ dbeta(1,0.5) in your model specification, and canopy <- frag$canopy and jags.params = "median" in your r code.

Answer 2

I think you could use a logit model for your probabilities. Maybe something like the following.

First, I convert your canopy observations back to the format that I suspect they began in, i.e. the number of canopy hits out of 20 samples at each site. I set 0.0001 to 0 and 0.999 to 1, and multiply the other canopy values by 20.

d$hits <- ifelse(d$canopy < 0.05, 0, ifelse(d$canopy > 0.95, 20, d$canopy * 20))

M <- function() {
  for (i in 1:n) {
    hits[i] ~ dbin(p[site[i]], 20)
  }
  for (j in 1:nsites) {
    logit.p[j] ~ dnorm(mu, sigma^-2)
    logit(p[j]) <- logit.p[j]
  }
  mu ~ dnorm(0, 0.0001) # uninformative prior for grand mean of logit(p)
  sigma ~ dunif(0, 10) # uninformative prior for sd of logit(p)
}

j <- jags(list(site=d$site, hits=d$hits, n=nrow(d), nsites=length(unique(d$site))), 
          NULL, 'p', M)

plot(j$BUGSoutput$summary[-1, '50%'], pch=20, xlab='site', xaxt='n', las=1, 
     ylim=c(0, 1), ylab=expression("p (median" %+-% "95% credible interval)"))
segments(1:13, j$BUGSoutput$summary[-1, '2.5%'], 
         y1=j$BUGSoutput$summary[-1, '97.5%'])
axis(1, 1:13, 1:13)