qpbayes - comparing bayes factor

[biopitch]

Dear Evan,

I am running qpbayes now on 15 qpbrute graphs. It seem to have run to completion but I get an empty heatmap pdf file. I am wondering if the script have ran properly and how to interpret the results.

I have attached my logs. qpbayes.log

Any advice is greatly appreciated.

[ekirving]

From that log file, it appears that everything ran to completion, so I don't know why the heatmap is empty. It could be an issue with one of the R packages used for plotting, but that's just a guess. Does the PDF raise an error when you open it? Is the filesize empty?

The heatmap is only intended as a visual representation of the Bayes factor analysis, and the raw data is available in the bayes/test1-likelihoods.csv file. It's also recommended to inspect each of the MCMC outputs, but most importantly the Gelman-Rubin metric, which should not exceed 1.2 for any of the models.

[biopitch]

Thanks for getting back to me so quickly.

The pdf opens with a notice that says 'The document contains no pages', the file size is 3.6kB.

Regarding the csv, I have something that looks like this: "","mean","sd" "e4497f6d8ac1",-266.332116919005,1.91874192054904 . . . "ff40788164f7",-298.824032288417,3.09625781106755

Does that mean 'ff40788164f7' is the best fit graph? But when I look at the Gelman-Rubin metric, many of the models are way above 1.2. Does that mean the results are not conclusive?

test-ff40788164f7-burn-gelman.pdf

[ekirving]

From the CSV, the model with the mean likelihood closest to 0 (i.e. the topmost graph) is the best supported.

However, there are two major caveats to consider:

1. The Bayes factors are only valid if the Gelman-Rubin metrics are below 1.2. The plot you attached shows the relationship between the Gelman-Rubin metric and a longer potential burn-in of the model (the point estimates of this metric are at the end of the log file for each model, and this is what you should use to evaluate each model). What your PDF indicates is that for that particular model, the MCMC chain length and burn-in was insufficient to reach the standing distribution of the model for some parameters. The best way to deal with this is by increasing these values using the --iterations and --burnin flags. If you keep the same output folder then qpBayes will resume where it left off from the existing chains, to save having to run them again. To make an appropriate choice, look at all the *-burn-gelman.pdf plots to decide how far you need to burn-in the existing chains (to get the Gelman-Rubin metrics low enough), and then extend the total chain length by the same quantity.

2. Once all your models have suitable Gelman-Rubin metrics, you can then use the CSV to determine which was the best model. However, you cannot reject one model in favour of another unless the difference between their log-likelihood scores is above a certain threshold (see https://en.wikipedia.org/wiki/Bayes_factor#Interpretation)

[biopitch]

Many thanks for the detailed explanation. Please let me just clarify some things:

1. Does this means all the plots that starts with 'a' and 'edge' should have a 'shrink factor' below 1.2? Let's say all my *-burn-gelman.pdf are similar to this one, does that mean I should increase the burnin to around 600000 where the 'shrink factor' drops off drastically?
2. Once I have obtained suitable metrics, and I start to look at the CSV file. For example, the first two entries are: "","mean","sd" "e4497f6d8ac1",-266.332116919005,1.91874192054904 "8b0ebb187bae",-266.340680648362,1.87853073386484 The difference between the two is so low, does that neither models can be rejected? So only models with differences that are >10 can we reject the other model.

[ekirving]

1. Under ideal conditions, it is best when all parameters have a 'shrink factor' below 1.2. However, this is not always obtainable in practice, and the most important parameter to consider is the likelihood itself. For some models, the point estimate for one or more admixture proportions does not have a substantial impact on the likelihood of the model itself, which can cause the MCMC sampler to fail to converge for that specific parameter (resulting in an elevated Gelman-Rubin metric). In many cases, increasing the chain length and burn-in can fix this, but with some difficult models the ideal condition is unobtainable. I recommend you try increasing the chain length and burn-in first, but if this does not work, then a principled way of resolving the problem is to calculate the multivariate version of the Gelman-Rubin metric. This gives you a single summary statistic for all parameters (instead of one for each). If this multivariate statistic is below 1.2 then you are okay to proceed with the next step.
2. Given the negligible difference between the two models you showed above (K=0.0086) you have no power to reject the second model in favour of the first. I would generally use a minimum K value of ~3, but the larger this value is, the more certainty you have about rejecting the alternative model. However, until your models have converged, the mean likelihood values can be quite misleading, and I have seen situations where a model with a large Gelman-Rubin metric ultimately leapfrogged all the other models once it had properly converged. So I would advice extreme caution about interpreting the likelihood scores from models that have not converged.

[biopitch]

Thanks Evan for the clarification.

I am currently trying to run:

qpBayes \
--geno test.geno \
--ind test.ind \
--snp test.snp \
--prefix test \
--pops M St B J K \
--out Out \
--iterations 500000 \
--burnin 100000 \
--threads 12

The gelman-rubin graphs didn't change/get replaced, resulting in the same likelihood calculations. Do I first delete the existing files in the bayes folder?

[ekirving]

Both your --iterations and --burnin parameters are well below the defaults, in light of which, I am unsurprised that your models have not fully converged. Running a model to completion with the default parameters does take quite a while, however, the defaults were chosen specifically to avoid situations like this.

Producing the burn-gelman.pdf plots (and other diagnostics) also takes a long time, so by defualt, they are not regenerated if the files already exist. To remake them, delete all the output in the bayes/ folder except the chains files themselves.

[biopitch]

Thanks Evan, I will try that.

[ekirving]

FYI, I just release a new version of qpBayes with a few improvements to the way the Bayes factors are calculated (and a bug fix for the empty heatmap issue that you reported).

To see the benefit you will need reinstall qpBrute, and then delete the contents of both the dstats/ and bayes/ folders and run qpBayes again.

pip install --upgrade --no-deps --force-reinstall git+git://github.com/ekirving/qpbrute.git

The new version of the {prefix}-likelihoods.csv file contains the multivariate version of the Gelman-Rubin metric (i.e., the MPSRF) and {prefix}-bayes-factors.csv contains the pairwide K values for the Bayes factors.