[R-sig-ME] Prediction, Poisson and Intercepts
Manuel Morales
Manuel.A.Morales at williams.edu
Fri Oct 23 18:37:21 CEST 2009
Dear list,
I have a question about how to efficiently generate predictions from a
poisson mixed model in lmer. A previous post by Bert Gunter
(http://markmail.org/message/22ncrgbt76bj5bmo) suggested the following:
model.matrix(terms(a.lmer.model),new.data) %*% fixef(a.lmer.model)
However, the back-transformed values from these predictions do not match
the original means because the fixed effect estimate of the intercept
does not include the random effects (see code below).
One alternative is:
new.Intercept <- log(mean(exp(coef(a.lmer.model)[[1]][,1])))
new.fixed <- fixef(a.lmer.model)
new.fixed[1] <- new.Intercept
model.matrix(terms(a.lmer.model),newdata) %*% new.fixed
But this seems awfully clunky ... Any suggestions for a more elegant
solution?
Demonstration code for mean model below:
library(reshape)
library(lme4)
seed=1
data <- melt(sapply(c(.1,.5,1,5),function(x) rpois(50,x)))
mean(data$value) #Overall mean
m.glm <- glm(value~1,data=data, family=poisson)
exp(coef(m.glm)) #Exponent of intercept = population mean
m.lmer <- lmer(value~(1|X2),data=data, family=poisson, nAGQ=20)
ranef(m.lmer)$X2+fixef(m.lmer)==coef(m.lmer)
mean(exp(unlist(coef(m.lmer)))) #Mean of exp of coefs = overall mean
--
http://mutualism.williams.edu
More information about the R-sig-mixed-models
mailing list