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

[Mikhail Matz]

Hi Jarrod - 

I see how the means of random effects don't have to be normal (they just happen to be close to normal in my case, which is rather surprising). What I am concerned about is the Poisson residuals, corresponding to the difference between exp(latent variable) (latent=as.vector(apply(model$Liab,2,mean))) and the observed counts.  Each count should be a sample from a Poisson distribution with the mean=exp(corresponding latent variable), correct? This is what I still cannot convince myself about. Or is the sample truncated in some way (say, only allowing the numbers larger than or equal to the mean?)  

cheers

Misha
 

On Aug 29, 2012, at 10:09 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:

> Hi,
> 
> The residuals (pre observation random effects) are assumed to be conditionaly normal. However, I think that even if this is satisfied then this does not imply that the posterior means/modes of the random effects will be normal. In fact when the expected number of counts is small I could imagine that the posterior means/modes could be strongly non-normal, perhaps even multimodal. Are the "latents" the posterior mean/mode of the residuals?  If you plot the distribution of per-observation  random-effects on an iteration by iteration basis, do you still see non-normality?
> 
> 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.
> 
>