[R-sig-ME] lmer, poisson family and offset variable
Ivar Herfindal
ivar.herfindal at bio.ntnu.no
Wed Feb 7 17:24:36 CET 2007
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).
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"
More information about the R-sig-mixed-models
mailing list