VMCMC: a graphical and statistical analysis tool for Markov chain Monte Carlo traces
© The Author(s) 2017
Received: 7 November 2015
Accepted: 28 January 2017
Published: 10 February 2017
MCMC-based methods are important for Bayesian inference of phylogeny and related parameters. Although being computationally expensive, MCMC yields estimates of posterior distributions that are useful for estimating parameter values and are easy to use in subsequent analysis. There are, however, sometimes practical difficulties with MCMC, relating to convergence assessment and determining burn-in, especially in large-scale analyses. Currently, multiple software are required to perform, e.g., convergence, mixing and interactive exploration of both continuous and tree parameters.
We have written a software called VMCMC to simplify post-processing of MCMC traces with, for example, automatic burn-in estimation. VMCMC can also be used both as a GUI-based application, supporting interactive exploration, and as a command-line tool suitable for automated pipelines.
VMCMC is a free software available under the New BSD License. Executable jar files, tutorial manual and source code can be downloaded from https://bitbucket.org/rhali/visualmcmc/.
KeywordsConvergence Markov chain Monte Carlo Metropolis-Hastings Phylogenetics Software Visualization
Bayesian inference using Markov chain Monte Carlo (MCMC) is today a common and trusted approach in molecular phylogenetics (see, e.g., in [1–4]), and is used by phylogeny inference software such as MrBayes , BEAST , BAli-Phy , PrIMe  and JPrIME . A typical goal in phylogenetics is to determine evolutionary relationships, for a set of species or for genes of interest, but researchers may also be interested in other parameters, for example parameters related to substitution patterns or duplication/loss processes. An advantage with Bayesian phylogenetic inference is that you can obtain posterior distributions of evolutionary parameters, conditional on your data, where the evolutionary parameters can be classified as discrete parameters (e.g., phylogenetic trees) or continuous (e.g., duplication and loss rates) . For phylogenetics, Bayesian inference is often complicated by the continuous parameters’ dependency on discrete parameters such as trees, that due to their structure can have problems with mixing.
The output from MCMC, the trace, is a conditional sampling of unknown parameters, needing post-processing to yield the desired posterior distributions. Users need to inspect the trace for possible non-convergence, estimate burn-in (how many of the first samples should be discarded?), assess mixing (does it look like a random sampling of the posterior?), and compute parameter statistics. Post-processing can, however, be constrained by practical difficulties, in particular due to presence of both discrete and continuous parameters in the trace, and due to the need to manually extract useful information from MCMC runs, currently using multiple software (for example Tracer , AWTY  and CODA ). User-friendly interfaces supporting graphical and interactive parameter exploration currently has potential to improve.
When converged, MCMC samples are drawn from the underlying stationary distribution and the trace represents a sample from the posterior distribution of analyzed parameters. After removing burn-in samples, a converged chain displays little correlation between remaining samples, which also indicates good mixing. Convergence is, in theory, guaranteed in MCMC if the chain is allowed enough iterations. However in practice, it is not possible to estimate how many iterations are needed, nor is it possible to determine whether the chain has converged. Hence, heuristics are used to assess non-convergence and the initial “burn-in”. Software like CODA , AWTY , and Tracer  have support for this, and more, and for example MrBayes  has a convergence diagnostic based on standard deviation of split frequencies built in. These heuristics analyze a single parameter of the chain at a time and the decision of convergence is, usually, left to visual analysis of the trace and the heuristic value proposed for the parameter. For example, Tracer  suggests the heuristic that the effective sample size (ESS)  is greater than 200, for all parameters, as well as manual analysis of the trace to ascertain that the chain has indeed converged. Please note that the recommendation/heuristic of ESS greater than 200 by Tracer does not have a theoretical justification or a systematic study to support it and using a sufficiently large ESS for assessing convergence for high-dimensional MCMC is recommended by Gong and Flegal . In the convergence analysis of large data sets, e.g., for genome-wide analysis, the manual inspection of convergence becomes a bottleneck and one would like to rely on automated non-convergence assessment and burn-in estimation based on multiple convergence diagnostics. Therefore, we see a need for simplified large-scale convergence analysis.
There are several reasons to scrutinize individual traces. First, one may doubt automatic convergence/burn-in assessment, e.g., when different convergence diagnostics give contradictory assessments. Second, there may be doubts regarding the MCMC proposals, in particular, when the sampled posterior distribution is multimodal and seemingly converged multiple runs provide very different samples. Third, there may be questions regarding mixing of an MCMC chain, e.g., when a trace is stuck in a local optimum for a long time and the proposal acceptance ratios of parameters are low. Fourth, surprising results, e.g., when a new result contradicts previous studies, a trace should be scrutinized to establish that the MCMC has not behaved strangely. Mixing of parameters is an important aspect of a successful MCMC run, in particular for phylogeny analysis where one wants to make sure that the tree posterior is sampled well. Although trace statistics can give important indications, various visualizations are often needed to better understand the different aspects of a chain and its trace. We have experimented with visualizations that focus on trees and their relation to the trace, that we have not found in Tracer, AWTY, and CODA, and believe the community can benefit from them.
We present VMCMC, a tool for phylogenetic MCMC analysis, with support for analysis and exploration of chain convergence, burn-in estimation, trace visualization, parameter estimation, graphical display of parameter traces, which can run both as a command-line tool and as an application with a graphical user-interface (GUI).
VMCMC presents MCMC traces statistically and visually, and enables both automatic and interactive analysis. VMCMC also supports output from many popular MCMC programs, e.g., PrIMe, JPrIME, BEAST, and MrBayes.
Large-scale automated analysis: Applying MCMC on genome-wide data, with perhaps hundreds or thousands of phylogenies to estimate, users typically make arbitrary burn-in decisions and, e.g., throw away the initial 25% of a trace. VMCMC can be used on the command line, easily integrated in scripts, and can be asked to estimate a burn-in, assess non-convergence, and estimate parameters automatically.
Detailed trace exploration: When trying to understand a single MCMC trace, and decide whether it is plausible that it has converged, it can be illuminating to see where in a trace certain trees are found, determine majority rule consensus trees and other parameters for different parts of a trace, and have tree space visualized. VMCMC supports a high degree of interactivity for such exploration.
Simultaneous handling of both real and tree parameters: Currently, there is no software that simultaneously analyzes both tree and continuous parameters of a chain. The available MCMC analysis software either handle the continuous parameter (CODA and Tracer) or the tree parameter (AWTY) but not both simultaneously. VMCMC supports analysis of both types of parameters simultaneously and can, e.g., show the occurrence of a specific tree topology in the trace of a real parameter or show the effect of changing burn-in on both the real and tree parameter simultaneously.
Graphical trace visualization: VMCMC employs a user-friendly GUI to display the trace of a selected real parameter, which helps in seeing trends for a particular parameter as well as visually inspect whether this parameter may have converged (Fig. 1).
Tree analysis: VMCMC presents a posterior distribution of trees, sorted by probability and hence easy access to the tree with maximum a posteriori probability and the closest alternatives. Selected trees are shown in Newick format and are visualized (Fig. 2) in a panel using code from forester . A majority rule consensus tree can be computed, also restricted to selected trees, thus simplifying comparative tree analysis.
Tree space analysis: VMCMC calculates unweighted Robinson-Foulds distance  between each pair of trees and displays a two-dimensional projection (Fig. 3) using a multi-dimensional scaling technique [17–19], where similar trees are expected to cluster. Point sizes proportional to probabilities inform the user of the estimated posterior.
Visualizing tree mixing and parameter dependence: VMCMC can indicate where a selected tree is found in the trace, thus help visualizing tree mixing and how parameters and trees are correlated (Fig. 4).
Convergence and burn-in: Convergence assessment and burn-in estimation can be performed using Geweke’s convergence diagnostic , Gelman-Rubin’s convergence diagnostic , and Höhna-Sahlin’s ESS-based estimator [22, 23].
Analyzing parallel chains: Currently, traces from two parallel chains can be analyzed and visualized together (Fig. 5). Convergence can be tested on numerical posteriors by applying a Mann-Whitney U test, and tree split distributions  can be compared using a chi-square test for two samples. VMCMC provides the flexibility to perform parallel chain analysis using user-specified burn-ins for both chains. The parallel chains are appended to each other after removing specified burnins for both chains, and statistics and posterior distribution is calculated for joint traces. This is a useful feature in cases with poor convergence or high computational demands.
Command line tool: VMCMC can be used both as a GUI application and as a command line tool suitable as a component when writing scripts. Command line output is given in JSON format for easy parsing.
VMCMC tries to hide the computational costs inherent in the problem and utilizes multithreading to perform computations (e.g., in loading of graphics, and in calculation of statistics and of convergence diagnostics) independently from the user interface. Some calculations are delayed: the “Tree Analysis” tab is not prepared until the user asks for it to avoid unnecessary calculation of the needed, and expensive, two-dimensional Robinson-Foulds distance matrix. Furthermore, in order to show trees with non-negligible posterior frequency only and to speed up the loading of analysis on “Tree Analysis” tab, we have limited the computation of distance matrix to unique trees with posterior probability at least 0.2%. While VMCMC can work with any number of unique trees in the posterior, a warning for expecting a delay in loading the tab is generated if the posterior contains more than 45 unique trees with frequency greater than 0.2%.
Execution time from specifying the input file to presentation of graphical results in VMCMC
Size of tree
Results and discussion
For an example application of VMCMC, we used a tetrapod 18S ribosomal RNA dataset, named “DS1”, collected by Hedges et al. . DS1 is known to be problematic for convergence and mixing of the tree parameter [26, 27]. We ran JPrIME twice with the same settings (using the JC69 model and default parameters; the species tree was dated using TimeTree.org  and is available as Additional file 1), each time for 10 million iterations and sampling every 100 iterations, yielding two traces containing 105 samples which we name trace 1 and trace 2 (available in Additional files 2 and 3). Here we present how VMCMC can be used to evaluate mixing of the tree parameter for trace file (Additional file 3). The tutorial for VMCMC is available as Additional file 4.
Figure 1 shows the trace plot of DLRModelDensity parameter (this is the probability reported from the DLR submodel ) for trace 1, where the fuzzy caterpillar-like trace indicates that the trace does not show non-stationary behavior for DLRModelDensity parameter (i.e., the chain does not seem to be “stuck” in a state, nor still continuously improving). In the “Tree” tab, we see a summary of the tree posterior on the left. As an example, we have selected the top five (by frequency) trees and computed a majority-rule consensus-tree for them (Fig. 2). The consensus tree has five edges with support values lower than 1.0, indicating where the selected trees differ. Looking at the “Tree analysis” tab (Fig. 3), the size of each topology (represented by a colored circle) is proportional to the frequency of each topology in the trace. The “Tree analysis” tab also shows that the MCMC chain is repeatedly sampling T0 and similar trees (e.g., T2 and T5), or T1 and similar trees (e.g., T3 and T4), but the transition to other regions is rare. The clear partitioning of tree space into two disjoint parts is indicative of a bimodal distribution. This conclusion is in agreement with previous studies [26, 27] that suggest the data is bimodal. If we select the five top trees, push the “Mark selection in graph” button, and view the “Graph” tab again (Fig. 4), the regions where the selected trees are found in the trace are marked in different colors. We observe intermixing of all five colored lines and no long undisturbed stretches of a single color.
Comparing two parallel chains, run independently on the same data, is another way to investigate MCMC mixing. VMCMC can load two traces and superimpose plots. Figure 5 displays superimposed trace plots, from the DS1 dataset, where one can see that the red trace is different from blue trace at many places. Such aberrations are also checked through statistical hypothesis testing, and we note that the Mann-Whitney U test (left sidebar in Fig. 5) rejects the hypothesis that the traces are from the same distribution and a two-sample chi-square test on trees’ split frequencies indicates the same (not shown here). We conclude that the two parallel chains are sampling from two different posterior distributions and have either not reached the stationary distribution, or need more iterations to be able to safely conclude convergence.
VMCMC has implemented various convergence diagnostics commonly used by phylogenetics community, e.g., Gelman-Rubin’s convergence diagnostic , effective sample size (ESS)-based diagnostics and Geweke convergence diagnostic . Please note that some of these diagnostics are outdated and these will either be removed or be replaced in the newer releases of VMCMC by more accurate and recent diagnostics like Fixed-Width Stopping Rule (FWSR) . To elucidate this point, Gong and Flegal  have shown that Geweke convergence diagnostic  is misleading and outdated for assessing convergence. Further calculations based on Gong and Flegal’s work  reveal that an ESS of 6000 measured by a 95% confidence interval corresponds to computational uncertainty of approximately 5% the size of the posterior standard deviation while the ESS recommendation of 200 made by Tracer and measured by a 95% confidence interval is equivalent to having a computational uncertainty of 28% the size of the posterior standard deviation. Therefore, the estimated variances used in calculations with Tracer’s heuristic can be unstable for the sample size of 200 and the simulation has, probably, not converged yet. Thus, we are looking forward to connect ESS based tests available in VMCMC to computational uncertainty using FWSR in the newer releases of VMCMC.
VMCMC can be applied to trace files from several molecular phylogenetics MCMC tools. Assessing whether a chain has converged and is sampling from the stationary distribution is a non-trivial task. As demonstrated by the example, VMCMC can help identify issues with mixing and convergence of the MCMC run for all parameters. The graphical user interface supports interactive data exploration and the command line interface enables large-scale automated application.
VMCMC simplifies tasks in MCMC analysis that we encounter in our work and we believe that our software can be valuable to the community as well.
Effective sample size
Fixed-width stopping rule
Graphical user interface
Markov chain Monte Carlo
The authors thank Kristoffer Sahlin for help on convergence statistics and burn-in estimation methods.
RHA was funded by a stipend from the higher Education Commission (HEC) of Pakistan. SAM was funded by EuroSPIN (an Erasmus Mundus joint doctoral program) and Swedish e-Science Research Center (SeRC). These funding bodies did not influence our design study, analysis, data interpretation, or manuscript.
Availability of data and materials
An executable jar file, tutorial, and source code of VMCMC can be downloaded from https://bitbucket.org/rhali/visualmcmc/. The datasets analysed during this study are included in this published article and its supplementary information files.
The first prototype was developed as Bachelors thesis project by JM and MB guided by LA. Further development has been done by RHA as part of his PhD project and he has been responsible for design, implementation, drafting the manuscript and the tutorial, and coordination between other co-authors. SAM implemented the parallel chain analysis part. RMA implemented the handling of input generated from MrBayes and BEAST software. JS has been involved with the tree analysis tab and Robinson-Foulds distance metric. SMZ has incorporated tree visualization code from forester. LA has supervised the project, generated the idea and drafted the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Ali RH, Khan AA. Tracing the evolution of FERM domain of Kindlins. Mol Phylogenet Evol. 2014; 80:193–204.View ArticlePubMedGoogle Scholar
- Albertin CB, Simakov O, Mitros T, Wang ZY, Pungor JR, Edsinger-Gonzales E, Brenner S, Ragsdale CW, Rokhsar DS. The octopus genome and the evolution of cephalopod neural and morphological novelties. Nature. 2015; 524(7564):220–4.View ArticlePubMedPubMed CentralGoogle Scholar
- Eyun S-i, Moriyama H, Hoffmann FG, Moriyama EN. Molecular evolution and functional divergence of trace amine–associated receptors. PloS ONE. 2016; 11(3):1–24.View ArticleGoogle Scholar
- Safavi-Hemami H, Li Q, Jackson RL, Song AS, Boomsma W, Bandyopadhyay PK, Gruber CW, Purcell AW, Yandell M, Olivera BM, Ellgaard L. Rapid expansion of the protein disulfide isomerase gene family facilitates the folding of venom peptides. Proc Natl Acad Sci. 2016; 113(12):3227–32.View ArticlePubMedPubMed CentralGoogle Scholar
- Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Hohna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012; 61(3):539–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007; 7:214.View ArticlePubMedPubMed CentralGoogle Scholar
- Suchard MA, Redelings BD. BAli-Phy: simultaneous Bayesian inference of alignment and phylogeny. Bioinformatics. 2006; 22(16):2047–8.View ArticlePubMedGoogle Scholar
- Åkerborg O, Sennblad B, Arvestad L, Lagergren J. Simultaneous Bayesian gene tree reconstruction and reconciliation analysis. Proc Natl Acad Sci USA. 2009; 106(14):5714–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Sjöstrand J, Sennblad B, Arvestad L, Lagergren J. DLRS: gene tree evolution in light of a species tree. Bioinformatics. 2012; 28(22):2994–5.View ArticlePubMedGoogle Scholar
- Ali RH. Phylogenetic inference. Chapter in Doctoral thesis From Genomes to Post-processing of Bayesian inference of Phylogeny, KTH Royal Institute of Technology Stockholm. 2016.
- Nylander JA, Wilgenbusch JC, Warren DL, Swofford DL. AWTY (are we there yet?): a system for graphical exploration of MCMC convergence in Bayesian phylogenetics. Bioinformatics. 2008; 24(4):581–3.View ArticlePubMedGoogle Scholar
- Plummer M, Best N, Cowles K, Vines K. CODA: Convergence diagnosis and output analysis for MCMC. R News. 2006; 6(1):7–11.Google Scholar
- Thiébaux HJ, Zwiers FW. The interpretation and estimation of effective sample size. J Clim Appl Meteorol. 1984; 23(5):800–11.View ArticleGoogle Scholar
- Gong L, Flegal JM. A practical sequential stopping rule for high-dimensional Markov Chain Monte Carlo. J Comput Graph Stat. 2016; 25(3):684–700.View ArticleGoogle Scholar
- Zmasek C. forester. https://github.com/cmzmasek/forester. Accessed 30 May 2016.
- Robinson DF, Foulds LR. Comparison of phylogenetic trees. Math Biosci. 1981; 53(1):131–47.View ArticleGoogle Scholar
- Kruskal JB. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika. 1964; 29(1):1–27.View ArticleGoogle Scholar
- Kruskal JB. Nonmetric multidimensional scaling: a numerical method. Psychometrika. 1964; 29(2):115–29.View ArticleGoogle Scholar
- Hillis DM, Heath TA, St John K. Analysis and visualization of tree space. Syst Biol. 2005; 54(3):471–82.View ArticlePubMedGoogle Scholar
- Geweke J. Bayesian Statistics 4. Oxford, UK: Clarendon Press: 1992.
- Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992; 7:457–511.View ArticleGoogle Scholar
- Höhna S. Burnin estimation and convergence assessment. Chapter in Licentiate thesis Bayesian Phylogenetic Inference, Stockholm University Stockholm. 2011.
- Sahlin K. Estimating convergence of Markov chain Monte Carlo simulations. Master’s thesis, Stockholm University. 2011.
- Huelsenbeck JP, Ronquist F. Bayesian analysis of molecular evolution using MrBayes In: Nielsen R, editor. Statistical Methods in Molecular Evolution. 1st. New York: Springer: 2005. p. 183–226.Google Scholar
- Hedges SB, Moberg KD, Maxson LR. Tetrapod phylogeny inferred from 18S and 28S ribosomal RNA sequences and a review of the evidence for amniote relationships. Mol Biol Evol. 1990; 7(6):607–33.PubMedGoogle Scholar
- Lakner C, Van Der Mark P, Huelsenbeck JP, Larget B, Ronquist F. Efficiency of Markov chain Monte Carlo tree proposals in Bayesian phylogenetics. Syst Biol. 2008; 57(1):86–103.View ArticlePubMedGoogle Scholar
- Whidden C, Matsen FA. Quantifying MCMC exploration of phylogenetic tree space. Syst Biol. 2015;006.
- Hedges SB, Marin J, Suleski M, Paymer M, Kumar S. Tree of life reveals clock-like speciation and diversification. Mol Biol Evol. 2015; 32(4):835–45.View ArticlePubMedPubMed CentralGoogle Scholar