ZIP MCMCglmm - how to increase effective sample size?
Hi Jarrod,
Thank you very much for your message!
I think one of the main issues is the high degree of autocorrelation. I tried to check for autocorrelation as Pierre suggested, but my computer ran out of memory.
Here are my answers to your questions:
1). There are 10523 observations with a mean of 0.09; 95% of observations are zero.
value N raw % valid % cumulative %
0 9946 94.52 94.52 94.52
1 372 3.54 3.54 98.05
2 112 1.06 1.06 99.12
3 58 0.55 0.55 99.67
4 22 0.21 0.21 99.88
5 7 0.07 0.07 99.94
6 4 0.04 0.04 99.98
7 1 0.01 0.01 99.99
8 1 0.01 0.01 100.00
missing 0 0.00
total N=10523 ? valid N=10523 ? x?=0.09 ? ?=0.44
2) group 1: 1768 level
group 2: 1857 levels
3) I attached the histogram of the latent variable.
4) I would like a model with two random intercepts (group 1 and 2) and I am not sure how to do that with a zero inflation term.
Thank you so much for taking the time to help me with this.
Best regards,
DNM
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> on behalf of Jarrod Hadfield <j.hadfield at ed.ac.uk>
Sent: Wednesday, November 1, 2017 1:42 PM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] ZIP MCMCglmm - how to increase effective sample size?
Sent: Wednesday, November 1, 2017 1:42 PM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] ZIP MCMCglmm - how to increase effective sample size?
Hi, I would try and diagnose why the mixing is bad first. 1) How many observations are there, what is their mean and how many are zero? 2) How many levels are there to group1 and group2? 3) What is the range of the latent variable for the zero inflation? Include pl=TRUE in the call to MCMCglmm and plot a histogram of model$Liab[,(ncol(model$Liab)/2+1):ncol(model$Liab)]) 4) why do you have random effects for zero-inflation but choose to ignore it in the fixed effects? Cheers, Jarrod On 01/11/2017 20:31, dani wrote: > Thanks, Pierre! Will do that, too! > > Best, > DNM > > ________________________________ > From: Pierre de Villemereuil <pierre.de.villemereuil at mailoo.org> > Sent: Wednesday, November 1, 2017 1:25 PM > To: dani > Cc: r-sig-mixed-models at r-project.org > Subject: Re: [R-sig-ME] ZIP MCMCglmm - how to increase effective sample size? > > You can look at the auto-correlation to guess how many more iterations are needed to increase your sample size. Something like: > > library(coda) > autocorr.diag(mj$Sol) > autocorr.diag(mj$VCV) > > Cheers, > Pierre > > On Wednesday, 1 November 2017 20:14:15 NZDT dani wrote: >> Hello Pierre and list members, >> >> >> Thank you so much! The analysis with the new prior worked:) However, the effective samples are still small, so I am trying again the new prior with more iterations - will report back how my effective samples change. >> >> >> Best regards, everyone! >> >> DNM >> >> >> ________________________________ >> From: Pierre de Villemereuil <pierre.de.villemereuil at mailoo.org> >> Sent: Tuesday, October 31, 2017 10:11 PM >> To: dani >> Cc: r-sig-mixed-models at r-project.org >> Subject: Re: [R-sig-ME] ZIP MCMCglmm - how to increase effective sample size? >> >> Just some parentheses issue: >> priori <- list(R=list(V=diag(2), nu=1, fix=2), >> G=list(G1=list(V=diag(2)/2, nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000), >> G2=list(V=diag(2)/2, nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000))) >> >> Cheers, >> Pierre >> >> Le mercredi 1 novembre 2017, 17:18:05 NZDT dani a ?crit : >>> Hi Pierre, >>> >>> I tried using the new prior you suggested and I got this error: >>> >>> Error in MCMCglmm(y ~ trait - 1 + at.level(trait,1):(x1+: >>> prior list should contain elements R, G, and/or B only >>> >>> I am not sure what to do about this:) >>> Any advice would be very much appreciated. >>> >>> Thanks, >>> DaniNM >>> <http://aka.ms/weboutlook> >>> >>> >>> ________________________________ >>> From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> on behalf of Pierre de Villemereuil <pierre.de.villemereuil at mailoo.org> >>> Sent: Tuesday, October 31, 2017 2:08 PM >>> To: r-sig-mixed-models at r-project.org >>> Subject: Re: [R-sig-ME] ZIP MCMCglmm - how to increase effective sample size? >>> >>> Hi, >>> >>> There are about three way to increase effective sample size: >>> - increase the number of iterations >>> - use a prior with better properties >>> - change your model somehow (you might not always want to use that one...) >>> >>> In your case, using a slightly more informative prior and the extended parameters prior might help? Something like: >>> >>> priori <- list(R=list(V=diag(2), nu=1, fix=2), >>> G=list(G1=list(V=diag(2)/2, nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000)), >>> G2=list(V=diag(2)/2, nu=2, alpha.mu=c(0,0), alpha.V=diag(2)*1000)))) >>> >>> Hope this helps, >>> Pierre. >>> >>> On Wednesday, 25 October 2017 18:45:46 NZDT dani wrote: >>>> Dear list members, >>>> >>>> I need some advice regarding this ZIP MCMCglmm model: >>>> >>>> library(MCMCglmm) >>>> >>>> priori <- list(R=list(V=diag(2), nu=0.002,fix=2), >>>> G=list(G1=list(V=diag(2), n=2),G2=list(V=diag(2), n=2))) >>>> >>>> mj <- MCMCglmm(y ~ trait - 1 + at.level(trait,1):(x1+x2+x3+x4+ x5 +x6+x7+ offset), >>>> random = ~idh(trait):group1 + idh(trait):group2, >>>> family = "zipoisson", >>>> prior = priori, >>>> rcov = ~idh(trait):units, >>>> verbose=FALSE, >>>> thin = 100, >>>> burnin = 3000, >>>> nitt = 103000, >>>> saveX=TRUE, saveZ=TRUE, saveXL=TRUE, pr=TRUE, pl=FALSE, >>>> data = s25h) >>>> >>>> summary(mj) >>>> >>>> # Iterations = 3001:102901 >>>> # Thinning interval = 100 >>>> # Sample size = 1000 >>>> # >>>> # DIC: 4811.791 >>>> # >>>> # G-structure: ~idh(trait):group1 >>>> # >>>> # post.mean l-95% CI u-95% CI eff.samp >>>> # traity.group1 0.4307 0.1351 0.9281 10.17 >>>> # traitzi_y. group1 4.3196 2.1216 7.4310 31.26 >>>> # >>>> # ~idh(trait):group2 >>>> # >>>> # post.mean l-95% CI u-95% CI eff.samp >>>> # traity. group2 0.4233 0.2341 0.6781 30.81 >>>> # traitzi_y. group2 3.5497 1.2365 6.1525 26.39 >>>> # >>>> # R-structure: ~idh(trait):units >>>> # >>>> # post.mean l-95% CI u-95% CI eff.samp >>>> # traity.units 0.02393 0.002833 0.06621 10.58 >>>> # traitzi_y.units 1.00000 1.000000 1.00000 0.00 >>>> # >>>> # Location effects: y ~ trait - 1 + at.level(trait, 1):(x1 + x2 + x3 + x4 + x5 + x6 + x7 + offset) >>>> # >>>> # post.mean l-95% CI u-95% CI eff.samp pMCMC >>>> # traity -4.3823820 -6.1496186 -2.6424402 23.592 <0.001 *** >>>> # traitzi_y 3.4696204 2.6430476 4.1392235 1.922 <0.001 *** >>>> # at.level(trait, 1):x1 -0.0498043 -0.2192051 0.1097667 16.979 0.522 >>>> # at.level(trait, 1):x2M -0.2088408 -0.4535085 0.0440055 8.727 0.088 . >>>> # at.level(trait, 1):x31 0.1422342 -0.1473884 0.4199985 11.521 0.288 >>>> # at.level(trait, 1):x4 0.0007054 -0.0030953 0.0043456 24.299 0.680 >>>> # at.level(trait, 1):x5 0.1131704 0.0647184 0.1676469 26.621 <0.001 *** >>>> # at.level(trait, 1):x6 -0.0128734 -0.0483344 0.0306350 13.599 0.588 >>>> # at.level(trait, 1):x7 0.0102356 -0.0276141 0.0540893 40.746 0.680 >>>> # at.level(trait, 1):offset 1.3511873 0.6963525 2.1299075 13.216 <0.001 *** >>>> # --- >>>> # Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 >>>> >>>> I would like to increase my effective samples, but I am not sure which way to go. I tried increasing the NITT to 503000, but the effective samples actually got worse. Is there anything else I could do? I plan on dropping some variables from the model, but if I were to proceed with the model above, what could I have done better? >>>> >>>> Thanks in advance! >>>> DNM >>>> <http://aka.ms/weboutlook> >>>> >>>> [[alternative HTML version deleted]] >>>> >>> _______________________________________________ >>> R-sig-mixed-models at r-project.org mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> stat.ethz.ch Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... > R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> stat.ethz.ch Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... > stat.ethz.ch > Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... > > > >> R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> stat.ethz.ch Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... >> stat.ethz.ch >> Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... >> >> >> > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-mixed-models at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> stat.ethz.ch Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models R-sig-mixed-models Info Page - stat.ethz.ch<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models> stat.ethz.ch Your email address: Your name (optional): You may enter a privacy password below. This provides only mild security, but should prevent others from messing ... -------------- next part -------------- A non-text attachment was scrubbed... Name: HistogtamLambda.png Type: image/png Size: 3205 bytes Desc: HistogtamLambda.png URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20171102/23ce076b/attachment-0001.png>