[R-sig-ME] lmer, poisson family and offset variable

[Douglas Bates]

On 2/7/07, Ivar Herfindal <ivar.herfindal at bio.ntnu.no> wrote:
> Dear Mixed-Models-list

> In poisson models, it can often be convenient with an offset variable,
> in order to account for known linear relationships that influence the
> counts (number of occurrences), e.g. observation time or patch size. In
> glm, (and, I hope, lmer), this is simply added as e.g.
> offset(log(exposure.time)) in the formula. This seems to work fine for
> me, and the parameter estimates, se and z-statistics seems reasonable.
> However, when I try the mcmcsamp to look at the distribution of the
> estimates, it looks horrible. Adding exposure time as fixed factor (i.e.
> log(exposure.time)) in the model instead of as an offset gives
> distributions that are reasonable and as expected. Do anyone know if
> offset is compatible with lmer and/or mcmcsamp, or can explain why the
> distribution become so weird when using exposure time as offset variable
> compared to as a fixed variable? Am I doing something completely wrong
> here?

> I provide a short example below (the data is perhaps not very realistic
> in terms of values, but provide a a design close to the one I am using.
> It works as an example to show you how the distribution looks like with
> exposure time as offset or fixed variable).

Thank you for providing the example.

I am not surprised that there is a problem with such an example.  I
can think of two possible causes - either I forgot to take into
account the offset when formulating the mcmcsamp or this could be
related to a generic problem with mcmcsamp for generalized linear
mixed models in which it gets stuck at particular values of the fixed
effects.  The mcmcsamp function applied to linear mixed models uses
exact sampling for three subsets of the parameters (\sigma^2, the
relative variance-covariance matrix for the random effects, and the
fixed effects and random effects sampled jointly).  I have used the
same approach in mcmcsamp applied to a generalized linear mixed model
but there is no exact sampling for the last group of parameters.  I
use a Metropolis-Hastings  method for those and sometimes it is
unsuccessful in moving the values of these parameters.

I would suggest looking at the time series plot of the MCMC sample
(see an earlier response on this list today for doing this with
library(coda); xyplot(...))  to see if the problem is the series
getting stuck.

> Otherwise, the lme4 package works wonderful and is the package I use
> most frequent in R. Thanks a lot to the developer(s).
>
> Sincerely
>
> Ivar Herfindal
>
> ##
> set.seed(100)
> testdata <- cbind.data.frame(
>     Obs=rpois(50, 20),
>     gr=rep(1:10, 5),
>     fixed=rnorm(50))
> testdata$exp.time <- testdata$Obs - (rnorm(50, 1,1))
> mod.offset <- lmer(Obs ~ offset(log(exp.time)) + fixed + (1|gr),
> data=testdata, family=poisson)
> mod.offset.mcmcsamp <- mcmcsamp(mod.offset, 10000)
> mod.no.offset <- lmer(Obs ~ (log(exp.time)) + fixed + (1|gr),
> data=testdata, family=poisson)
> mod.no.offset.mcmcsamp <- mcmcsamp(mod.no.offset, 10000)
>
> # histogram of parameter distribution with offset
> windows()
> par(mfrow=c(2,2))
> for(i in 1:3)
> hist(mod.offset.mcmcsamp[,i], col="gray",
> main=colnames(mod.offset.mcmcsamp)[i])
>
> # histogram of parameter distribution without offset
> windows()
> par(mfrow=c(2,2))
> for(i in 1:4)
> hist(mod.no.offset.mcmcsamp[,i], col="gray",
> main=colnames(mod.no.offset.mcmcsamp)[i])
>
> ##
> sessionInfo()
> R version 2.4.1 (2006-12-18)
> i386-pc-mingw32
> attached base packages:
> [1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
> "methods"   "base"
>
> other attached packages:
>        lme4      Matrix     lattice
> "0.9975-11"  "0.9975-8"   "0.14-16"
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>