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,
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])
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
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).
I think that the glmmADMB package is exactely what you are looking for.
install.packages("glmmADMB", repos="http://r-forge.r-project.org")
But in a bayesian point of view you can use the MCMCglmm package or the BUGS/JAGS software, they are very flexible and you can fit this kind of model. (and the syntax is close to the R one)
EDIT thanks to @randel
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")
install.packages("R2admb")
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.
