下面我使用负二项式，因为它比简单的泊松模型更灵活。数据是计数$y$16 个人的活动$x$. 每个人有 14 个计数（即计数周期）。

似然函数是dnegbin($p_x , r_x)$用一个单独的$p$和$r$对于每个人$x$从伽马分布中提取。（实际上我已经根据Kruschke重新参数化了这个位。）我遵循 Kruschke 对 gamma 先验的选择$r$和$m$这里。

现在我的问题是如何为负二项式的参数来自的两个伽马分布的参数选择超先验：

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

我已经尝试过伽玛分布dgamma(a,b)和设置$a$和$b$0.5 或 1.0 似乎都能产生合理的结果。

1. 为伽马分布的两个参数选择伽马先验是否合理，即两者都应该r.a来自r.b伽马分布？
2. 在这种情况下，信息最少的伽马先验的参数是什么，例如dgamma(0.5,0.5)？

.

JAGS/BUGS 型号：

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
}
"

.

运行模型

为了实验，我提供了一些数据：

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)

运行 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)

为每个人提取参数（随机效应）：

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)

绘制结果：

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())

对于任何正在努力将 JAGS 和 R 的负二项式论点联系起来的人：

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