Hyper-prior for negative binomial in hierarchical model using JAGS/BUGS

Question

Below I'm using a negative binomial because it is more flexible than a simple poisson model. The data are counts $y$ of events for 16 individuals $x$. There are 14 counts (i.e. counting periods) for each individual.

The likelihood function is dnegbin($p_x , r_x)$ with a separate $p$ and $r$ for each individual $x$ drawn from a gamma distribution. (Actually I have re-parameterized this bit according to Kruschke.) I follow Kruschke's choice of gamma priors for $r$ and $m$ here.

Now my question is how to pick the hyper-priors for the parameters of the two gamma distributions from which the parameters of the negative binomial come:

r[j] ~ dgamma(r.a,r.b)
m[j] ~ dgamma(m.a,m.b)

I've experimented with gamma distributions dgamma(a,b) and setting $a$ and $b$ both to 0.5 or 1.0 seems to yield reasonable results.

1. Is the choice of gamma priors for both parameters of a gamma distribution reasonable, i.e. should both r.a and r.b come from a gamma distribution?
2. What are the parameters of a minimally informative gamma prior in this case, e.g. dgamma(0.5,0.5)?

.

The JAGS/BUGS model:

modelString = "
model {
    for( i in 1:length(y) ) {
      y[i] ~ dnegbin( p[x[i]] , r[x[i]] )
    }
    for( j in 1:Nj ) {
        p[j] <- r[j]/(r[j]+m[j])
        r[j] ~ dgamma(r.a,r.b)
        m[j] ~ dgamma(m.a,m.b)
        v[j] <- r[j]*(1-p[j])/(p[j]*p[j])
    }
    r.a ~ dgamma(a, b)
    r.b ~ dgamma(a, b)
    m.a ~ dgamma(a, b)
    m.b ~ dgamma(a, b)
    a <- 1.0
    b <- 1.0
}
"

.

Running the model

For experimentation I provide some data:

d <-
list(x = rep(1:16, each=14),
y = c(3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 4, 0,
2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 2, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
3, 0, 0, 0, 0, 0, 4, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 2, 0,
0, 3, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 0, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), Nj = 16)

Run JAGS:

require(rjags)
parameters <- c("m","r","p","v","r.a","r.b","m.a","m.b")
jagsModel <- jags.model(textConnection(modelString), data=d, n.chains=3, n.adapt=1e3)
update(jagsModel, n.iter=1e3)
codaSamples <- coda.samples(jagsModel, variable.names=parameters, n.iter=3000, thin=1)
m <- as.matrix(codaSamples)

Extract parameters (random effects) for each individual:

out <- data.frame()
for (p in c('m','r','p','v')) {
    for (i in 1:16) {
        out[i,'x'] <- i
        c <- sprintf('%s[%i]',p,i)
        #out[i,p] <- mean(estimates[c])        # mean
        dens <- density(m[,c])
        out[i,p] <- dens$x[which.max(dens$y)]  # mode
    }
}

# calculate frequencies for each individual 'x' and #occurrences 
t <- data.frame(table(d$x,cut(d$y, breaks=0:6, right=F)))
colnames(t) <- c('x','breaks','freq')

# add fits to frequency table
t <- merge(t,out)
t <- within(t, freq.fit <- dnbinom(as.numeric(breaks)-1, size=r, mu=m)*14)

Plot results:

library(ggplot2)
ggplot(t,aes(breaks,freq,group=x)) + geom_point() + geom_line(aes(y=freq.fit)) + facet_wrap(~x) + theme(panel.grid=element_blank(), panel.background=element_blank())

For anyone who is struggling to relate the arguments of the negative binomials from JAGS and R:

JAGS,  R's dnbinom
   r,         size
   p,         prob
   m,           mu