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

[Manuel Morales]

On Fri, 2009-10-23 at 18:04 +0100, Jarrod Hadfield wrote:
> 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.

Trying that in the demo code below doesn't seem to work for, but maybe
I'm misunderstanding (or more likely I'm not being clear).

Using the variance of the single random effect:
exp(fixef(m.lmer)+.5* 2.3184) equals 2.093279. 

When I try mean(exp(unlist(coef(m.lmer)))) or mean(data$value), both
equal 1.595.

> 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
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> 
> 
-- 
http://mutualism.williams.edu