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

[Jarrod Hadfield]

Hi Manuel,

Sorry - I misunderstood.  I thought you wanted to predict what the  
mean would be for some fixed effects, under a new random set of random  
effects - rather than the 4-levels of X2 you observed. The function I  
sent does the former.

Cheers,

Jarrod



Quoting Manuel Morales <Manuel.A.Morales at williams.edu>:

> 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
>
>
>



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