xavierdidelot / BactDating

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

Sigma value 501.000 on multiple runs on same dataset! #41

Closed maesaar closed 3 years ago

maesaar commented 3 years ago

I was wondering if there is maximum sigma value? I have started multiple runs on the same dataset: 1) Run with correct dates 2) Forced dates 3) Permuted dates runs 1 - 6 (to check the clock rate)

I have got different results for mu, sigma and alpha. But sigma=501.000000 occurres multiple times. For example in the first run with correct dates and in multiple runs with permuted dates.

I have two questions: 1) Is the sigma value 501.000000 just coincidence or is something wrong with my results? 2) Is it reasonable to compare models (DIC value) only when both models has converged (all three ESS values > 200) - because sometimes I have got results when model with all three converged ESS values have bigger DIC than model which some of ESS values have not converged (same as in #40)?

Thank you!

maesaar commented 3 years ago

Additional clarification: the forced dates model did not converge (with same parameters used as original model). In this case is it correct to compare DICs of two models?

maesaar commented 3 years ago

One more question I have ran three separate runs on same data with same parameters: 1) Correct dates: res=bactdate(t,d,nbIts=2000000,model="mixedcarc",useRec=T,showProgress=T) 2) Forced dates: res2=bactdate(t,dforce,nbIts=2000000,model="mixedcarc",useRec=T,showProgress=T)

First run results: res (mu=295.1867; sigma=501.0000; alpha=249.8216 and DIC=2430.68) res2 (mu=1.135386; sigma=395.243993; alpha=1.0808848 and DIC=4348.76)

Second run results: res (mu=206.8084; sigma=435.1820; alpha=288.6447 and DIC=2399.03) res2 (mu=3.361895; sigma=501.0000; alpha=8.673352 and DIC=2439.13)

Third run results: res (mu=215.3106; sigma=501.0000; alpha=183.9929 and DIC=2440.26) res2 (mu=6.069539; sigma=501.0000; alpha=5.065853 and DIC=2428.48)

Why is that the DIC is so different between runs?

dnanto commented 3 years ago
2. Is it reasonable to compare models (DIC value) only when both models has converged (all three ESS values > 200)

Yes, you shouldn't take models with low ESS values seriously as they exhibited poor Markov chain mixing. Also, see here: [https://mc-stan.org/docs/2_21/reference-manual/effective-sample-size-section.html]().

maesaar commented 3 years ago

@dnanto but does this mean that the model with not converged ESS is definitely not better despite of better DIC?

dnanto commented 3 years ago

@mmaesaar I don't understand the math of DIC vs ESS 100% yet. However, a model must converge to be taken seriously. This means that one should compare those with ESS > threshold.

Although, in my opinion, it's not as simple as picking the model with the best information criterion, whether it's DIC, AIC, BIC, etc... For instance, if a more complex model slightly outperforms a simpler one, then one should probably abide by Occam's Razor and go with the simpler one as long as ESS is satisfactory of course.

maesaar commented 3 years ago

Ok perhaps @xavierdidelot can enlighten us! Because I want to compare 1) correct dates model with forced dates to detect significant temporal signal; 2) correct dates model with permuted dates model for temporal signal.

Now is the question of ESS convergence vs DIC value

xavierdidelot commented 3 years ago

When the ESS is too small the results should not be used, even for model comparison. I have no idea why you get repeated 501 values for sigma, it's not a limitation in BactDating itself but might be inherited from other packages or a natural limit for the computation. In any case, with such a high sigma the clock is so relaxed that it becomes meaningless and all dating becomes more or less random. I think your problem stems from the fact that you are trying to run with randomised or forced-equal dates, so that the temporal signal is non-existent which makes the model very hard to fit, because the MCMC needs to explore a vast range of dates. With randomised dates you could look at the value of mu instead of DIC, and if you find that it is smaller with randomisation than without then this is good evidence for the significance of the temporal signal, cf

Duchene S, Duchene D, Holmes EC, Ho SY. 2015. The performance of the date-randomization test in phylogenetic analyses of time-structured virus data. Mol. Biol. Evol. 32:1895–1906.

maesaar commented 3 years ago

@xavierdidelot thank you for the quick explanation. Did I understand correctly that if the runs with correct dates sigma is 501.000 that means the dating results are meaningless? - So it is not advised to use results for example for the run:

res (mu=295.1867; sigma=501.0000; alpha=249.8216 and DIC=2430.68)

xavierdidelot commented 3 years ago

Yes, I'm afraid the results are probably meaningless in this case.

maesaar commented 3 years ago

I have many ranges in my dataset perhaps this is the cause - I will refine my dataset to specific dates then re-run analysis.

Lastly, if I compare only mu's (for example correct dates vs forced dates and correct dates vs permutated dates) then I should check that mu is smaller and CIs of mu between different runs don't overlap?

xavierdidelot commented 3 years ago

Yes that's right, you should check that the CI of permuted run is lower and does not overlap with the CI when using correct dates.

maesaar commented 3 years ago

Thank you!