Skip to content

MCMCglmm poisson / not poisson

4 messages · Mikhail Matz, Jarrod Hadfield

Original
Mikhail Matz
# 
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
Jarrod Hadfield
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi,

If you calculate latentsv (a vector of variances from the posterior  
distribution of the residuals), and replace

exp(latents[i])

with

exp(latents[i]+0.5*latentsv[i])

does this improve things?

Cheers,

Jarrod

Quoting Mikhail Matz <matz at utexas.edu> on Tue, 28 Aug 2012 22:06:45 -0500:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      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: