Prediction, Poisson and Intercepts

Hi Manuel,

The function below is taken from the MCMCglmm course notes that will  
be available soon. It  generates predictions for GLMM, and the random  
effects can be marginalised if you wish (I think this is what you  
want?). For the Poisson distribution this can be done exactly, for  
distributions such as the binomial, approximations have to be used.

mu is the fixed effect prediction, var is the sum of the variance  
components you want to average over.

Cheers,

Jarrod



predict.MCMCglmm<-function(mu, var=NULL, family="gaussian"){
   if(family=="gaussian"){
     return(mu)
   }
   if(family=="poisson"){
     return(exp(mu+0.5*var))
   }
   if(family=="binomial"){
     return(inv.logit(mu-0.5*var*tanh(mu*(1+2*exp(-0.5*var))/6)))
   }
}



On 23 Oct 2009, at 17:37, Manuel Morales wrote:

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