[R-sig-ME] ZIP MCMCglmm

Jarrod Hadfield j.hadfield at ed.ac.uk
Sun Oct 7 01:07:49 CEST 2012


Hi,

A) The residual variance for the zero-inflation is not identified nor  
is the residual covariation between the zero-inflation and the Poisson  
part. However, you are attempting to estimate both. Use:

rcov=idh(trait):units

and fix the [2,2] element of the covariance matrix at something (e.g  
1) in the prior (i.e. add fix=2 to your R prior).

B) you have no unique intercepts for the zero inflation and the  
Poisson part. You should have trait as a main effect in the fixed  
formula (or trait-1 might be easier to interpret).

C)

  random= ~us(1+sqkm):EFF_Track

assumes that the regression parameters (intercept and slope) of sqkm  
on log Poisson rate and logit zero-inflation are identical within  
EFF_Tracks, and that the slope=0 on average (because its not in the  
fixed formula). I would have one of the following in the fixed effects:

a) trait:sqkm
b) at.level(trait,1):sqkm
c) at.level(trait,2):sqkm

and one of the following in the random effects.

i) random= ~us(trait):EFF_Track
ii) random= ~idh(trait):EFF_Track
iii) random= ~us(at.level(trait,1)):EFF_Track
iv) random= ~us(at.level(trait,2)):EFF_Track

and make sure that if i) or ii) to have a),  if iii) to have a) or b)  
and if iv) to have a) or c).

If you use ZAP rather than ZIP then fixed=~sqkm and random= ~EFF_Track  
is OK too.

D) DIC should not be trusted.

E) Not sure, but it sounds like g(0) is correcting for missed  
sightings of individuals present (?) and therefore refers to the  
Poisson part of the model.  The zero-inflation part may also be due to  
missed sightings of individuals present (e.g. on some transects the  
observers had their eyes closed!) but it is possibly more likely to be  
due to absence of porpoises? I don't know.

Cheers,

Jarrod









Quoting Katrin Ronnenberg <ronnenberg at gmx.com> on Fri, 05 Oct 2012  
12:31:12 +0200:

> Dear all,
> I have several questions on a zero-inflated MCMCglmm model.
>
> Data consist of 764 transects with porpoise counts ("sum_ppho"). On  
> 304 transects zero porpoises were seen. Data come from several  
> studies from the North Sea, they were all sampled in roughly the  
> same area and from altogether 12 years. We are interested in a trend  
> ("year"), seasonal effects ("dayofyear" quadratic (centered) Term +  
> the simple term. Moreover we have a simple geographical variable  
> "we"; a factor with two levels "West" and "east". The transect ID  
> (EFF_Track) I integrated as random factor, assuming that transects  
> are not independant of each other I used the us structure. Most  
> transects were revisited over several years, but not all of them.  
> Since Transects had differing length and visibility conditions,  
> which affected the effective area sampled, I used the effective area  
> ("sqkm") as random slope term.
>
>
>
> in MCMCglmm I tried the following:
>
> prior<-list(R=list(V=diag(2),nu=0.002), G=list(G1=list(V=diag(2),nu=0.002)))
>
> Mz<-MCMCglmm(sum_ppho~at.level(trait,1):year*at.level(trait,1):we+
> +           at.level(trait,1):(I(cdayofyear^2)+  
> at.level(trait,1):dayofyear), random= ~us(1+sqkm):EFF_Track,
> +              family="zipoisson",
> +              data=encounter, rcov=~us(trait):units, prior=prior,  
> verbose=FALSE,
> +              nitt=100000,thin =1000, burnin = 5000)
>> summary(Mz)
>
> Iterations = 5001:99001
>  Thinning interval  = 1000
>  Sample size  = 95
>
>  DIC: 2843.102
>
>  G-structure:  ~us(1 + sqkm):EFF_Track
>
>                                   post.mean l-95% CI u-95% CI eff.samp
> (Intercept):(Intercept).EFF_Track   3.62368  1.21791  6.03721    6.215
> sqkm:(Intercept).EFF_Track         -0.38983 -0.61860 -0.15801   14.857
> (Intercept):sqkm.EFF_Track         -0.38983 -0.61860 -0.15801   14.857
> sqkm:sqkm.EFF_Track                 0.05025  0.02699  0.07369   42.211
>
>  R-structure:  ~us(trait):units
>
>                               post.mean  l-95% CI u-95% CI eff.samp
> sum_ppho:sum_ppho.units          0.8150  0.673652   1.0178   69.615
> zi_sum_ppho:sum_ppho.units       0.1894 -0.472509   0.6914    7.978
> sum_ppho:zi_sum_ppho.units       0.1894 -0.472509   0.6914    7.978
> zi_sum_ppho:zi_sum_ppho.units    0.1561  0.001532   0.6305    8.428
>
>  Location effects: sum_ppho ~ at.level(trait, 1):year *  
> at.level(trait, 1):we + at.level(trait, 1):(I(cdayofyear^2) +  
> at.level(trait, 1):dayofyear)
>
>                                     post.mean   l-95% CI   u-95% CI  
> eff.samp  pMCMC
> (Intercept)                        -3.315e+00 -4.486e+00 -2.045e+00   
>   1.993  <0.01 *
> at.level(trait, 1):year             1.673e-03  1.111e-03  2.169e-03   
>   3.976  <0.01 *
> at.level(trait, 1):wew             -5.459e+02 -6.696e+02 -4.612e+02   
>  40.273  <0.01 *
> at.level(trait, 1):I(cdayofyear^2) -3.117e-05 -5.063e-05 -7.107e-06   
>  95.000  <0.01 *
> at.level(trait, 1):dayofyear       -1.726e-03 -3.851e-03  2.335e-05   
>  59.665 0.0632 .
> at.level(trait, 1):year:wew         2.727e-01  2.305e-01  3.343e-01   
>  40.388  <0.01 *
> ---
>
>
>
>
> In general I feel insecure if my assumptions are right and whether  
> prior specifications make sense. But I have two main questions:
> 1. Question: Am I correct in the assumption, that, if the output for  
> autocorr(M1$Sol) and autocorr(M1$VCV) are low, that spatial  
> autocorrelation is no issue? If I specify a relatively simple model  
> (only main effects), autocorrelation is no issue, however the more  
> interaction terms I include the stronger the autocorrelation  
> becomes. Increasing itterations by a factor of ten to:  
> nitt=100000,thin =5000, burnin = 10000, hardly changed anything  
> about that. Is there anything else I could try? I know that I can  
> increase nitt... further and ask Malcolm about the sir() function.  
> But I'm especially interested, if experimenting with the random term  
> has a chance of improving things?
> 2. Question: Following methods of distance sampling, a G(0) factor  
> was estimated compensating for the diving porpoises, which escaped  
> the surveyors. Originally we corrected the observations by dividing  
> the counts by the G(0) factor. So, basically only the transects with  
> counts were corrected. Am I right in the assumption that the latent  
> variable in a ZIP model compensates for all values with false zeros  
> i.e. also the false zero for the third animal in a transect with 2  
> counted animals? If that's the case, I would just not include any  
> correction for g(0). But if the latent variable only compensates for  
> false zeros of all zero values, but not for other values from the  
> poisson distribution, then I would assume that the correction factor  
> is needed. Is that the difference between a hurdle and a ZIP model?  
> ZIP was in my case better than hurdle according to DIC.
>
> I'm sorry, I hope you can understand my questions. I'm neither a  
> statistician, nor a native speaker, so I would like to appologize  
> for bad choice of words. And please let me know, if a representive  
> sample of my data would help.
> Thank you for reading this, ideas and suggestions are very much appreciated.
>
> Katrin
>
> _______________________________________________
> 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