[R-sig-ME] Fwd: MCMCglmm poisson / not poisson

[Mikhail Matz]

I thought I better post the end of this discussion, for the benefit of future doubters…

Begin forwarded message:

> From: Mikhail Matz <matz at utexas.edu>
> Subject: Re: [R-sig-ME] MCMCglmm poisson / not poisson
> 
> Interesting!
> so each mcmc realization of latent variables is perfectly poisson-compatible, but their mean is not. OK! This totally works for me. 
> I now wonder why I expected the means to retain Poisson properties in the first place, actually
> 
> Thanks a lot!
> 
> Misha
> 
> On Aug 29, 2012, at 10:39 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
> 
>> Hi,
>> 
>> library(MCMCglmm)
>> l<-rnorm(100,1,2)
>> y<-rpois(100, exp(l))
>> 
>> dat<-data.frame(y=y)
>> 
>> m1<-MCMCglmm(y~1, data=dat, family="poisson", pl=TRUE)
>> 
>> pp.poisson(y, colMeans(m1$Liab))
>> 
>> # looks bad
>> 
>> pp.poisson(y, colMeans(m1$Liab)+0.5*apply(m1$Liab, 2, var))
>> 
>> # looks better
>> 
>> for(i in 1:1000){
>> pp.poisson(y, m1$Liab[i,])
>> }
>> 
>> # looks good
>> 
>> Cheers,
>> 
>> Jarrod
>> 
>> 
>> 
>> 
>> 
>> Quoting Mikhail Matz <matz at utexas.edu> on Tue, 28 Aug 2012 22:06:45 -0500:
>> 
>>> 
>>> Hello -
>>> 
>>> I am playing with ways to justify that the MCMCglmm model fits my data well, which is quite important for me since I am hoping to be able to suggest MCMCglmm-based modeling as a general solution for a particular type of analysis.
>>> 
>>> I am running "poisson" family on counts data, with two random effects. Following Elston, D. A., R. Moss, et al. (2001). Parasitology 122: 563-569., I am checking whether my lognormal residuals (latent variable minus predicted value) are normally distributed (check), if my random effects (saved with pr=T) are normally distributed (more or less check), and then I try to see if the observed counts really look like Poisson samples based on the latent variables. Again, following Elston et al, I am making a p-p plot using this script (expert coders, please don't judge):
>>> 
>>> pp.poisson=function(counts,latents) {
>>> 	sim=c()
>>> 	for(i in 1:length(counts)){
>>> 		if (is.na(counts[i])) next
>>> 		data=counts[i]
>>> 		low=ppois(data,exp(latents[i]))-dpois(data,exp(latents[i]))
>>> 		up=ppois(data,exp(latents[i]))
>>> 		ss=seq(low,up,(up-low)/100)
>>> 		sim=append(sim,sample(ss,1))
>>> 	}
>>> 	sims=sort(sim)
>>> 	xx=(rank(sims)-0.5)/length(sims)
>>> 	plot(sims~xx)
>>> 	abline(0,1)
>>> }
>>> 
>>> … and unfortunately it looks really ugly, like a very strongly bent  ' ~ ' rather than a line.
>>> The little script above seems to work; here is a sanity check:
>>> 
>>> psim=c()
>>> nnn=rnorm(500,10,10)
>>> for (i in 1:length(nnn)){
>>> 	psim=append(psim,rpois(1,exp(nnn[i])))
>>> }
>>> pp.poisson(psim,nnn)
>>> 
>>> I will be extremely grateful for any comments on this.
>>> 
>>> cheers
>>> 
>>> Misha
>>> UT Austin
>>> 
>>> _______________________________________________
>>> 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.
>> 
>> 
>