[R-sig-ME] ZIP MCMCglmm - how to increase effective sample size?
dani
orchidn at live.com
Wed Nov 1 05:18:05 CET 2017
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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list