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

[Jarrod Hadfield]

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


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.