[R-sig-ME] Prediction, Poisson and Intercepts

[Luca Borger]

Not sure this helps at all, but you get the same result for the mean by this:


> mean(fitted(m.lmer))
[1] 1.674732
> 


HTH


Cheers,

Luca


----- Messaggio originale -----
Da: "Manuel Morales" <Manuel.A.Morales at williams.edu>
A: "r-sig-mixed-models" <r-sig-mixed-models at r-project.org>
Inviato: Venerdì, 23 ottobre 2009 12:37:21 GMT -05:00 U.S.A./Canada, stati orientali
Oggetto: [R-sig-ME] Prediction, Poisson and Intercepts

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

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models