您可以使用 R 以多种方式将多级 GLMM 与泊松分布（过度分散）拟合。很少R有包是：lme4、、、MCMCglmm等arm。一个很好的参考是Gelman 和 Hill (2007)

我将给出一个使用rjagspackage in执行此操作的示例R。它是R和JAGS（如OpenBUGS或WinBUGS）之间的接口。

$$n_{ij} \sim \mathrm{Poisson}(\theta_{ij})$$ $$\log \theta_{ij} = \beta_0 + \beta_1 \mbox{ } \mathtt{Treatment}_{i} + \delta_{ij}$$ $$\delta_{ij} \sim N(0, \sigma^2_{\epsilon})$$ $$i=1 \ldots I, \quad j = 1\ldots J$$ $\mathtt{Treatment}_i = 0 \mbox{ or } 1, \ldots, J-1 \mbox{ if the } i^{th} \mbox{ observation belongs to treatment group } 1 \mbox{, or, } 2, \ldots, J$

上面代码中的部分模拟了过度分散。但是没有人阻止你对个体之间的相关性（你不相信个体是真正独立的）和个体内部的相关性（重复测量）进行建模。此外，速率参数可以通过一些其他常数进行缩放，如. 请参阅 Gelman 和 Hill (2007) 以获取更多参考。这是简单模型的代码：$\delta_{ij}$rate modelsJAGS

data{
        for (i in 1:I){         
            ncount[i,1] <- obsTrt1[i]
            ncount[i,2] <- obsTrt2[i]
                ## notice I have only 2 treatments and I individuals 
    }                               
}

model{
    for (i in 1:I){ 
        nCount[i, 1] ~ dpois( means[i, 1] )
        nCount[i, 2] ~ dpois( means[i, 2] )

        log( means[i, 1] ) <- mu + b * trt1[i] + disp[i, 1]
        log( means[i, 2] ) <- mu + b * trt2[i] + disp[i, 2]

        disp[i, 1] ~ dnorm( 0, tau)
        disp[i, 2] ~ dnorm( 0, tau)

    }

    mu  ~ dnorm( 0, 0.001)
    b   ~ dnorm(0, 0.001)
    tau ~ dgamma( 0.001, 0.001)
}

这是R实现使用它的代码（比如说它被命名为overdisp.bug：）

dataFixedEffect <- list("I"       = 10,
                        "obsTrt1" = obsTrt1 , #vector of n_i1
                        "obsTrt2" = obsTrt2,  #vector of n_i2
                        "trt1"    = trt1,     #vector of 0
                        "trt2"    = trt2,     #vector of 1
                       )

initFixedEffect <- list(mu = 0.0 , b = 0.0, tau = 0.01)

simFixedEffect <- jags.model(file     = "overdisp.bug",
                             data     = dataFixedEffect,
                             inits    = initFixedEffect,
                             n.chains = 4,
                             n.adapt  = 1000)

sampleFixedEffect <- coda.samples(model          = simFixedEffect,
                                  variable.names = c("mu", "b", "means"),
                                  n.iter         = 1000)

meansTrt1 <- as.matrix(sampleFixedEffect[ , 2:11])
meansTrt2 <- as.matrix(sampleFixedEffect[ , 12:21])

你可以玩弄你的参数的后验，你可以引入更多的参数来让你的建模更精确（我们喜欢这样想）。基本上，你明白了。

有关使用rjagsand的更多详细信息JAGS，请参阅John Myles White 的页面

无需离开 lme4 包即可解决过度分散问题；只包括观察数的随机效应。提到的 BUGS/JAGS 解决方案可能对您来说太过分了，如果不是，您应该有易于拟合的 lme4 结果进行比较。

data$obs_effect<-1:nrow(data)
overdisp.fit<-lmer(y~1+obs_effect+x+(1|obs_effect)+(1+x|subject_id),data=data,family=poisson)

这在这里讨论：http : //article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 Elston 等人在非正式和学术上。（2001 年）。

我认为 glmmADMB 包正是您正在寻找的。

install.packages("glmmADMB", repos="http://r-forge.r-project.org")

但是从贝叶斯的角度来看，您可以使用MCMCglmm包或BUGS/JAGS软件，它们非常灵活，您可以适应这种模型。（并且语法接近 R 的）

编辑感谢@randel

如果你想安装glmmADMB和R2admb包，最好这样做：

install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)   
install.packages("R2admb")

到目前为止的好建议。这里还有一个。rhierNegbinRw您可以使用包的功能拟合分层负二项式回归模型bayesm。