[R-sig-ME] Question re: MCMCglmm with zipoisson family

[Jarrod Hadfield]

Dear Charlotte,

The structure of the model looks valid, although you may be able to  
simplify it more and get the same level of fit.

Initially, I would work with a reparameterisation of your model:

abund3.mcmc <-
MCMCglmm(abund~trait+trait:fl+trait:mps+trait:elev+trait:distance 
+trait:dens-1,
random=~idh(trait):YEAR 
+idh(trait):pop,rcov=~idh(trait):units,data=flyden
s,family="zipoisson",prior=priorD)

This fits the same model, but the interpretation is probably easier in  
this case.

The first fixed term will be the intercept of the Poisson process and  
the second fixed term will be the intercept of the zero-inflation. The  
remaining terms are trait (Poisson and Zero-inflation) specific  
contrasts from the intercept.

You may also want to fix the zero inflation as a constant across all  
levels (both in the fixed and random terms). For example,

priorB <- list(R=list(V=diag(2),n=2,fix=2),
G=list(G1=list(V=diag(c(1, 1e-6)),n=2, fix=2),G2=list(V=diag(c(1,  
1e-6)),n=2, fix=2)))

fixes the random effect variances of the zero-inflation to zero, which  
may be OK.  However, having 1 as your prior variance for the Poisson  
process should not be used by default, and you should make sure your  
conclusions don't depend on it.

Also, you may not want zero-inflation to vary across your fixed effect  
terms (e.g. fl).  All your fixed terms look continuous, so your  
parameters (using the above parameterisation) should be of the form:

1) poisson intercept
2) zi intercept
3) effect of fl on poisson process
4) effect of fl on zi process
5) effect of elev on poisson process
6) effect of elev on zi process

and so on...

you could set effects 4,6,8.... to zero by fixing them  at zero in the  
prior specification:


priorC <- list(B=list(mu=matrix(0, nterms, 1),  V=diag(nterms)*1e+6),  
R=list(V=diag(2),n=2,fix=2),
G=list(G1=list(V=diag(c(1, 1e-6)),n=2, fix=2),G2=list(V=diag(c(1,  
1e-6)),n=2, fix=2)))

# this sets up a diffuse prior around zero on the fixed effects, where  
nterms is the number of fixed effects your fitting

diag(priorC$B$V)[seq(4, nterms, 2)]<-1e-6

# this sets the variance around the zero inflation terms (except the  
intercept) to be very small, essentially fixing them to zero.


One potential problem you may have with diffuse priors on logit-binary  
outcomes is when the probability is extreme. In this  case the logit  
probability can be very large or very small and huge changes in the  
logit probability actually translate into very small differences on  
the probability scale. For example,

 > inv.logit(5)
[1] 0.9933071
 > inv.logit(10)
[1] 0.9999546

In this case it may make more sense to put an informative prior on the  
zi fixed effects. I often use a variance of pi^2/3 because this is  
approximately uniform on the probability scale:

hist(inv.logit(rnorm(10000, 0, sqrt(pi^2/3))))

but there are probably better options when there are random effects.  
Andrew Gelman has written something on this.

Cheers,

Jarrod



On 4 Aug 2009, at 15:55, Klank Charlotte wrote:

> Dear list,
>
>
>
> I am currently working on a data set on the fly abundance in flowers  
> of
> Trollius europaeus in 20 populations, recorded over three years.
>
>
>
> I am interested in the of effect of several explanatory variables such
> as plant population size, mean plant size, plant density, elevation  
> and
> distance to next sampled population on the fly abundance in a
> population.
>
>
>
> As I think that my data is zero-inflated (it contains 87% zero  
> values),
> I would like to use a zero-inflated model, which also allows me to fit
> fixed effects as well as random effects to account for the populations
> being repeatedly measured over three years.
>
>
>
> I therefore tried to fit the MCMCglmm with a zipoisson family and am
> wondering if any of you could give me some advice on the coding, as  
> I am
> not really sure if I have done things correctly.
>
>
>
> I have checked the MCMCglmm documentation as well as previous posts
> regarding MCMCglmm with a zipoisson family, but still have some  
> problems
> with how to code and if my priors are ok.
>
>
>
> So I would greatly appreciate it if you could have a look at the code
> below and let me know if it is ok and/or point me at the right reading
> material.
>
>
>
> ---
>
> Data set:
>
> n= 16867
>
> counts done in 2006,2007 and 2008
>
> 9 populations were sampled in 2006, all 20 in 2007 & 2008
>
>
>
> response:
>
> abund -  number of flies per flower, ranging from 0 - 7
>
>
>
> variables:
>
> YEAR -  year of sampling (factor)
>
> pop - name of the sampled population (factor)
>
> fl - population size (as No. of flowers) (int)
>
> elev - elevation (int)
>
> mps - mean plant size (average No. of flowers per individual) (num)
>
> dens -  plant population density (num)
>
> distance - distance to closest sampled population (num)
>
>
>
>
>
> MODEL CODE:
>
>
>
> priorA <- list(R=list(V=diag(2),n=2,fix=2),
> G=list(G1=list(V=diag(2),n=2),G2=list(V=diag(2),n=2)))
>
>
>
> abund3.mcmc <-
> MCMCglmm(abund~trait:fl+trait:mps+trait:elev+trait:distance 
> +trait:dens,r
> andom=~idh(trait):YEAR 
> +idh(trait):pop,rcov=~idh(trait):units,data=flyden
> s,family="zipoisson",prior=priorD)
>
>
>
> ----
>
>
>
> Thanks already for the help :>
>
>
>
> Regards,
>
>
>
> Charlotte
>
>
>
>
>
>
>
>
>
>
> 	[[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.