[R-sig-ME] Back-transformation of Poisson model

[Mollie Brooks]

Hi R-sig-mixed,
I found the below discussion of back-transformations from the log scale and the same formula is on page 46 of the MCMCglmm course notes pdf (version June 20, 2015). I can’t find a copy of the cited book McCulloch CE, Searle SR (2001) or any other source for this information. I have one point that I’m still wondering about…should v, the sum of the variance components, include residual error in any model, i.e. not just residual error that is part of an observation-level random effect for over-dispersed models? In the following example, the formula, exp(Xb+0.5*v) seems to be giving me an estimate that is consistently higher than my population mean. So adding the residual variance will make it even worse. Is this normal?
Now for the R-specific part of the question... is there a version of predict.merMod that gives the the predictions on the response scale averaged over possible realizations of the random effects, i.e. marginalized with respect to the random effects? I can’t get any predictions close to the observed mean in the following example. I expected that a new level of the random effect would be given a marginalized prediction. Maybe this could be added to the documentation.
thanks,
Mollie

library(lme4)

# simulate Poisson data with random intercept
set.seed(11)
dat=data.frame(group=paste0('g', rep(1:10, each=100)))
x=rnorm(10, mean=5,sd=2) 
dat$y = rpois(nrow(dat), exp(x[dat$group]))
mean(dat$y)#observed mean

# fit model
m1=glmer(y ~ (1 | group), family = poisson, dat)

#estract estimates
str(predict(m1,  type="response", re.form=NA))
exp(fixef(m1)+ 0.5*VarCorr(m1)$group[1])#from MCMCglmm Coursenotes p.46
#not very close to the mean

nd=data.frame(group=paste0('g', -5:-1)) #add new levels
str(predict(m1, newdata=nd, type="response", re.form=NULL, allow.new.levels=TRUE))
str(predict(m1, newdata=nd, type="response", re.form=NA, allow.new.levels=TRUE))
str(predict(m1, newdata=nd, type="response", re.form=NA))
str(predict(m1, type="response", re.form=NA))
mean(exp(coef(m1)[[1]][,1]))#only number close to the observed mean

> 
> ----- Forwarded message from j.hadfield at ed.ac.uk <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> -----
>      Date: Tue, 12 Jan 2010 18:38:06 +0000
>      From: Jarrod Hadfield <j.hadfield at ed.ac.uk <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>>
> Reply-To: Jarrod Hadfield <j.hadfield at ed.ac.uk <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>>
>   Subject: Re: [R-sig-ME] Back-transformation of Poisson model
>        To: "Kardynal,Kevin [Yel]" <Kevin.Kardynal at EC.gc.ca <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>>
> 
> Dear Kevin,
> 
> You will have to exponentiate the output from predict.
> 
> It is worth noting there are three types of prediction (on the
> data-scale) you could make:
> 
> a) exp(Xb)
> b) exp(Xb+Zu)
> c) int exp(Xb)du assuming the u's are iid.
> 
> X - fixed effect design matrix  b - fixed effects
> Z - random effect design matrix  u - fixed effects
> 
> You have done the first, but this can be very different from c, which
> is what I expect you want: the predictions averaged over possible
> realisations of the random effects.
> 
> To obtain c use exp(Xb+0.5*v)
> 
> where v is the sum of the variance components. If you use additive
> over-dispersed models you want to include the "residual" variance in
> v, also.
> 
> Cheers,
> 
> Jarrod
> 
> Quoting "Kardynal,Kevin [Yel]" <Kevin.Kardynal at EC.gc.ca <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>>:
> 
> > Hello,
> >
> > I'd like to know if predicted values from  a Poisson mixed model require
> > back-transformation to get the 'real' predicted values or if this is
> > done automatically? I assume that the predicted values are already
> > back-transformed since my results are similar to the annual means.
> >
> > My model for estimating bird trends in lme4 is as such:
> >
> > lme3<-lmer(Abundance ~ Year +  (1|Observer) + (1|Site/Station),
> > data=SWTH, family = poisson(link=log))
> >
> > I then took the coefficients from the mixed model and used them in a GLM
> > to derive the predicted values.
> > dd <-fixef(lme3) # gives coefficients
> >
> > lme_glmDEY <- glm(Abundance ~ Year, data=SWTH, family = poisson(link =
> > log)) # for poisson distribution
> > lme_glmDEY$coefficients <- dd
> >
> > # create a data set of predicted values
> > pred <- as.data.frame(predict(lme_glmDEY, SWTH , "response", se.fit=T))
> > pred_val <- as.data.frame(cbind(SWTH,pred))
> >
> > So, do I have to perform a back-transformation on the predicted values
> > (ie., using exp())?
> >
> > Thanks,
> >
> > Kevin
> >
> > 	[[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models <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.
> 
> 
> 
> ----- End forwarded message -----
> 
> 
> -- 
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
> 
------------------------
Mollie Brooks, PhD
Postdoctoral Researcher, Population Ecology Research Group
Department of Evolutionary Biology & Environmental Studies, University of Zürich
http://www.popecol.org/team/mollie-brooks/ <http://www.popecol.org/team/mollie-brooks/>
	[[alternative HTML version deleted]]