search

xavierdidelot / BactDating

Bayesian inference of ancestral dates on bacterial phylogenetic trees
https://xavierdidelot.github.io/BactDating
MIT License
80 stars 15 forks source link

ESS values when nbIts increase #53

Closed snakano0527 closed 2 years ago

snakano0527 commented 2 years ago

Hi,

I use BactDating for the first time. I ran bactdating with relaxedmodel for the same dataset (Gubbins output of 100 pneumococcal genome data) with two different nblts values, 1e5 and 1e7. I expected that all mcmc output values would increase when I ran with nblts=1e7 than ones with =1e5. However, the results were opposite; all values were lower when nblts=1e7 than when nblts=1e5. Is this possible? Or any setting errors?

please give me any advice.

Best regards, Satoshi

xavierdidelot commented 2 years ago

Hi Satoshi,

Longer runs will not necessarily result in larger ESS. In any case the ESS is capped at the number of sampled iterations, which is equal to nbIts/thin/2. Unless you specify a value for thin, this cap will be equal to 500 which is what happened in #41 (I think there was a confusion in this issue between the value of sigma and the ESS on sigma) It's always better to do longer runs if possible (so using nbIts=1e7 rather than nbIts=1e5), and if you get ESS values of at least 100 or 200 then the results can be trusted.

Best wishes, Xavier

amyecampbell commented 3 months ago

Hi @xavierdidelot ! Thank you for developing this tool. Sorry to reopen an old issue, but since my question is related, I thought I might ask it here. Why is it that the maximum effective sample size is 500 by default, and not 1000, given that the thinning interval is nbits/1000? In other words, where does the '2' come from in nbIts/thin/2?

amyecampbell commented 3 months ago

Hi @xavierdidelot ! Thank you for developing this tool. Sorry to reopen an old issue, but since my question is related, I thought I might ask it here. Why is it that the maximum effective sample size is 500 by default, and not 1000, given that the thinning interval is nbits/1000? In other words, where does the '2' come from in nbIts/thin/2?

AH. Is it due to the 50% burn-in?

xavierdidelot commented 3 months ago

Yes that's it, by default BactDating removes 50% as burn-in. You can change this in the function as.mcmc.resBactDating(x, burnin = 0.5)