[R-sig-ME] ZIP MCMCglmm - how to increase effective sample size?

[Pierre de Villemereuil]

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]]
>