[R-sig-ME] ZIP MCMCglmm - how to increase effective sample size?
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Nov 1 21:42:52 CET 2017
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>
>> 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
--
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