[R-sig-ME] Poor mixing and autocorrelation of ZIP data in MCMCglmm
Pierre de Villemereuil
p|erre@dev|||emereu|| @end|ng |rom ephe@p@|@eu
Thu Jun 11 18:17:23 CEST 2020
Hi,
Yes, you need to fix the "residual" variance of the ZI part to 1 and cannot estimate the "residual" covariance (as far as I know for the latter), but you still can estimate the variances of the random effects.
But this output doesn't seem right. The variance for "zi_baitIntake.id" (140!!) is way too high and the uncertainty around the fixed effects (especially the intercepts for each trait) is considerable.
Are you sure you have enough data to fit such a complicated model?
You can maybe be even more informative toward small values for the variances by setting nu = 1000 for the random effects part:
list(R = list(V = diag(2), nu = 2, fix = 2),
G = list(G1 = list(V = diag(2), nu = 1000, alpha.mu = c(0,0), alpha.V = diag(2)),
G2 = list(V = diag(2), nu = 1000, alpha.mu = c(0,0), alpha.V = diag(2)))
But PLEASE NOTE this prior above would not be valid if you wanted to estimate a covariance between the Poisson and ZI part for the random effects (i.e. if you used "us()" instead of "idh()" in your model definition).
The priors (previous and the above one) are not an inverse Gamma distribution, so do not state that. You can state the values for the parameters V, nu, alpha.mu and alpha.V, or even better, share the analysis code with your Methods.
Cheers,
Pierre
Le jeudi 11 juin 2020, 07:46:34 CEST Vital Heim a écrit :
> Dear Ben and Pierre,
>
> Thank you very much for the quick response - that is amazing and very
> helpful!
>
> I have not centered the continuous predictors as I had the feeling that
> this would not be necessary in my model. Would you recommend to center them?
>
> Thank you for the prior suggestion. From reading the course notes to
> MCMCglmm I thought that I have to fix all the random effects variances to 1
> in the ZI part. I am very new to Bayesian statistics and MCMCglmm and
> therefore used the MCMCglmm course notes to write and understand my model.
> There was one part in the course notes that said that there is no residual
> variance for the zero-inflated process and that we cannot estimate the
> residual covariance between the zero-inflation and the Poisson process but
> that we can deal with this by fixing the residual variance for the ZI
> portion at 1. Therefore I thought I have to do the same for my model.
>
> I used the new prior and the model ran fine and showed some improvement.
> The mixing of the random effects improved but the ZI portion of the random
> effects did not mix well. Do you think the mixing could be/needs to be
> further improved? Would you increase nu? I added you the model summary and
> the autocorrelation diagnostics below:
>
> > summary(model10)
>
> Iterations = 150001:1849501
> Thinning interval = 1500
> Sample size = 1134
>
> DIC: 4039.85
>
> G-structure: ~idh(trait):id
>
> post.mean l-95% CI u-95% CI eff.samp
> baitIntake.id 0.2097 0.01955 0.5446 1134.00
> zi_baitIntake.id 140.2144 10.39550 412.5714 16.88
>
> ~idh(trait):dive
>
> post.mean l-95% CI u-95% CI eff.samp
> baitIntake.event 0.2386 1.113e-01 0.3711 544.1
> zi_baitIntake.event 1.6603 1.054e-07 5.7216 64.6
>
> R-structure: ~idh(trait):units
>
> post.mean l-95% CI u-95% CI eff.samp
> baitIntake.units 0.5084 0.423 0.5999 1134
> zi_baitIntake.units 1.0000 1.000 1.0000 0
>
> Location effects: bait ~ trait - 1 + trait:sharks + trait:divers +
> trait:boats + at.level(trait, 1):presence + at.level(trait, 1):temperature
> + trait:gender
>
> post.mean l-95% CI u-95% CI eff.samp pMCMC
>
> baitIntake -0.989171 -8.506002 7.022529 1036.3 0.80071
>
> zi_baitIntake 5.726646 -0.731875 14.671429 118.8 0.09171
> .
> baitIntake:conspecifics -0.049999 -0.115503 0.015834 1134.0 0.13757
>
> zi_baitIntake:conspecifics 0.185816 -0.129883 0.601750 437.8 0.31922
>
> baitIntake:tourists 0.044321 0.010587 0.075382 1134.0 0.00529
> **
> zi_baitIntake:tourists -0.154735 -0.341183 0.021306 258.0 0.05996
> .
> baitIntake:operators 0.157988 0.069429 0.255920 1134.0 0.00176
> **
> zi_baitIntake:operators -0.232157 -0.776007 0.249766 472.4 0.35979
>
> presence 0.010091 0.008344 0.011620 1134.0 < 9e-04
> ***
> temperature 0.260196 -0.062494 0.536471 1027.0 0.08466
> .
> baitIntake:gendermale 0.355938 -0.329717 1.143065 1134.0 0.28571
>
> zi_baitIntake:gendermale -1.787142 -16.777974 10.576166 200.1 0.80952
>
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> > autocorr.diag(model10$Sol)
> baitIntake zi_baitIntake baitIntake:conspec zi_baitIntake:consp.
> baitIntake:tourists zi_baitIntake:tourists
> Lag 0 1.00000000 1.00000000 1.000000000 1.000000000
> 1.00000000 1.0000000000
> Lag 1500 0.04457168 0.32629918 0.011215824 0.400247111
> -0.02362611 0.4065199561
> Lag 7500 -0.01118381 0.21355704 -0.056422233 0.050903050
> -0.06575857 0.1543329508
> Lag 15000 -0.01182550 0.14366161 -0.004787545 0.048794655
> 0.01045932 0.1180500281
> Lag 75000 -0.06205008 0.05428778 -0.012689660 0.004470857
> -0.01587028 -0.0001966966
> baitIntake:operators zi_baitIntake:operators presence
> temperature baitIntake:gendermale zi_baitIntake:gendermale
> Lag 0 1.000000000 1.00000000 1.000000000
> 1.000000000 1.000000000 1.00000000
> Lag 1500 0.018702637 0.34990109 -0.005337675
> 0.049056906 -0.008540371 0.11797192
> Lag 7500 0.029576245 0.04236116 0.051737798
> -0.004486268 -0.043298182 0.10059317
> Lag 15000 -0.006467939 0.07262102 0.020471250
> -0.011604859 0.010096511 0.04657296
> Lag 75000 0.050869606 -0.01289837 0.034904725
> -0.064627241 -0.019433089 0.04871585
>
> > autocorr.diag(model10$VCV)
> baitIntake.id zi_baitIntake.id baitIntake.event
> zi_baitIntake.event baitIntake.units zi_baitIntake.units
> Lag 0 1.00000000 1.0000000 1.00000000
> 1.000000000 1.00000000 NaN
> Lag 1500 -0.01347780 0.6908849 0.03201474
> 0.789290686 0.02634913 NaN
> Lag 7500 -0.04617914 0.5013689 -0.03223097
> 0.519696189 -0.02148069 NaN
> Lag 15000 0.01414618 0.5099142 -0.04454962
> 0.341925156 -0.04412612 NaN
> Lag 75000 -0.03180043 0.2266740 -0.01645174
> -0.003025432 -0.01177512 NaN
>
> And then I would have another question regarding the prior which I am not
> sure if that is ok to post here. If not feel free to let me know. I have
> been looking up papers to learn how to properly report my prior choice and
> details in a publication. I found some papers that just wrote "we chose a
> inverse gamma prior" but I am concerned that some reviewers would want more
> specific information. Do you have a recommendation on how to correctly
> report a prior in a publication?
>
> Best wishes,
> Vital
> --
> *Vital Heim *
>
> *PhD student*
> University of Basel, Switzerland
>
> Bimini Biological Field Station Foundation
> Bimini, Bahamas
>
> +41 (0)79 732 05 57
> vital.heim using gmail.com
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list