# How do I fit a multilevel model for over-dispersed poisson outcomes?

I want to fit a multilevel GLMM with a Poisson distribution (with over-dispersion) using R. At the moment I am using lme4 but I noticed that recently the quasipoisson family was removed.

I've seen elsewhere that you can model additive over-dispersion for binomial distributions by adding a random intercept with one level per observation. Does this apply to the poisson distribution as well?

Is there a better way to do it? Are there other packages that you would recommend?

You can fit multilevel GLMM with a Poisson distribution (with over-dispersion) using R in multiple ways. Few R packages are: lme4, MCMCglmm, arm, etc. A good reference to see is Gelman and Hill (2007)

I will give an example of doing this using rjags package in R. It is an interface between R and JAGS (like OpenBUGS or 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$

The $\delta_{ij}$ part in the code above models overdispersion. But there is no one stopping you from modeling correlation between individuals (you don't believe that individuals are really independent) and within individuals (repeated measures). Also, the rate parameter may be scaled by some other constant as in rate models. Please see Gelman and Hill (2007) for more reference. Here is the JAGS code for the simple model:

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


Here is the R code to implement use it (say it is named: 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,

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


You can play around with your parameters' posteriors and you can introduce more parameters to make you modeling more precise (we like to think this). Basically, you get the idea.

For more details on using rjags and JAGS, please see John Myles White's page

• Thanks!! I only recently started looking into bayesian analysis and i still find it a bit difficult to grasp. I guess this is an opportunity to learn a bit more about it. Feb 9, 2011 at 9:16
• Why not gamma dispersion? Apr 18, 2011 at 9:14
• @Patrick you can definitely do that. But since I am taking log of the mean I prefer normal disp effect. A log normal distribution is another way to modeling distributions similar to gamma distribution. HTH. Apr 20, 2011 at 7:54

No need to leave the lme4 package to account for overdispersion; just include a random effect for observation number. The BUGS/JAGS solutions mentioned are probably overkill for you, and if they aren't, you should have the easy to fit lme4 results for comparison.

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)


This is discussed here: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 informally and academically by Elston et al. (2001).

• What if a model consists of two nominal variables, one continuous variable (all as fixed effects) and one grouping variable (random effect) with 3rd-order interactions and, moreover, number of measured subjects is equal to number of observations or records in the dataset? How should I cover this in the model? Jan 28, 2014 at 13:38

I think that the glmmADMB package is exactely what you are looking for.

If you want to install the glmmADMB and R2admb packages it is better to do:
install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)

• I believe currently the package should be installed via install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos") plus install.packages('R2admb'). Aug 4, 2014 at 15:53
Good suggestions so far. Here's one more. You can fit a hierarchical negative binomial regression model using the rhierNegbinRw function of the bayesm package.