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

Jarrod Hadfield j.hadfield at ed.ac.uk
Mon Mar 14 21:04:15 CET 2016 

Hi Mollie,

In simple models it is the sum of the variance components, although with 
things such as random regression it gets a bit more complicated.  For 
log-link Poisson models the mean is simply the mean of a log-normal. For 
logit models there is not an anlytical solution, and the approximation 
in the CourseNotes I certainly got from somewhere but can't now remember 
where from: somebody else has pointed out that my McCulloch & Searle 
citation does not have the relevant material so I will try and hunt down 
where the result comes from.

I think the deviation of your predicted and actual mean arises because 
of sampling error in the intercept/variance. The mean of an 
exponentiated sum of two normals will be larger than the exponentiated 
sum of the means of two normals. Consequently an unbiased estimator of 
the intercept and variance would lead to an upwardly biased estimator of 
the data-scale mean. Not sure about this though! Perhaps see if the 
predicted and actual mean become closer if you up the number of random 
effect levels (and so the intercept/variance become more precisely 
estimated).

Regarding Q2 I'm also not sure. In predict.MCMCglmm you can specifiy the 
random effects you want to marginalise over, but I *think* this is not 
implemented in lme4?

Cheers,

Jarrod







On 14/03/2016 16:34, Mollie Brooks wrote:
> 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]]
>
> _______________________________________________
> 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.

More information about the R-sig-mixed-models mailing list