jueves, 31 de marzo de 2016

Equivalence Between Maximum Parsimony and Maximum Likelihood

The maximum parsimony method is based on choosing the tree with the fewest changes character states (Camin Sokal, 1965). The maximum likelihood method is based on choosing the tree most likely on the data (observations) (Sober 2004). The debate over which method is best for inferring phylogenetic relationships has been in force. One of the arguments to justify the maximum parsimony method as valid, is that under some models is equivalent to the method of maximum likelihood (Farris, 1983). Tuffley and Steel (1997) showed that under non common mechanism model, select the same tree likelihood that selects parsimony. In this work, we compare the topologies selected by the parsimony and likelihood under a symmetrical model.

For this, use a tree of 10 tips. From this, generate sequences 20.000 bp in seq-gen (Rambaut Grassly, 1997) under the JC model and use 20 replicas. Trees reconstructed using parsimony method in TNT (Goloboff et al 2008) and in PHYML likelihood method (Guindon S et al, 2010) under a fully symmetrical model. Finally I contrast trees using the metric of Robinson-Foulds (Robinson et al 1981) as a measure of distance. This is done with the function treedist of phangorn package R.

The analysis got that for 20 replicas that do, only in a distance between the tree generated likelihood and generated parsimony is 0 (Fig 1, Table 1).

table 1. Symetric di erence for 20 replicas

Fig 1. Diagram of the distance between the trees generated by ML and MP

This is not the expected result, however, several parameters that I not considered, could influence in the result. One of them and I consider that most affected was the result was the type of method to use in the analysis of parsimony. Tuffley and Steel in his work used Fitch parsimony, I use Wag-
ner parsimony. This affects the order of the characters because the characters wagner considered as additives, while Fitch are considered as non-additives. Therefore, it could not conclude that the maximum likelihood method is not equivalent to the maxim parsimony method.


Camin J.H., Sokal R.R. 1965. A method for deducing branching sequences in a phylogeny. Evolution. 19:311–326

Sober, E. (2004). The contest between parsimony and likelihood. Systematic biology, 53(4), 644-653.

Farris, J. (1983). The logical basis of phylogenetic analysis.

Tuffley, C. and Steel, M. (1997): “Links Between Maximum Likelihood and Maximum Parsimony

under a Simple Model of Site Substitution.” Bulletin of Mathematical Biology 59: 581-607.

Rambaut A, Grassly NC, Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees., Comput Appl Biosci, June 1, 1997

Goloboff, P. A., Farris, J. S. and Nixon, K. C. (2008), TNT, a free program for phylogenetic analysis. Cladistics, 24: 774–786. doi: 10.1111/j. 1096-0031. 2008.00217.x

New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Guindon S., Dufayard J.F., Lefort V., Anisimova M., Hordijk., Gascuel O. 2010, Systematic Biology, 59(3):307-321

Robinson, D., Foulds, L. R. (1981). Comparison of phylogenetic trees. Mathematical bioscienes, 53(1), 131-147.

Dirichlet pior

                                                                        Sensibility analysis: Dirichlet prior
Compound Dirichlet priors for branch lengths was proposed as a solution to the problem of branch length overestimation in bayesian inference because the default exponential prior produces too large branch lengths and in some cases are unreal. Several authors have investigated this problem: Brown et al. (2010), Marshall (2010),  and Rannala et al. (2012).  Rannala et al. (2012) suggest the main cause is the poor default prior on branch lengths used in MrBayes, which assigns independient and identical distribution (i.i.d.) for branch lengths of the tree and proposed the compound Dirichlet priors on branch lengths, this uses a diffuse gamma or inverse gamma distribution on the tree length and the Dirichlet distribution is used to partition the tree length into branch lengths. The goal of this present work is compare default prior on branch lengths with the new compound Dirichlet priors implemented in MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003).

First, download sequences for 94 taxa of the froglet Crinia, genes : 12S, 16S, trnF, tRNA-Val using GenBank. Align by genes in MUSCLE (multiple sequence comparison by log-expectation) (Edgar, R.C. 2004). Set the evolutionary model using AIC in the software R version 3.2.3 (R Core Team, 2014) using phangorn (Schliep, 2011) package and the function modelTest. Create a nexus file to run  the Bayesian analysis in MrBayes version 3.1.2 (Ronquist and Huelsenbeck 2003), for compound Dirichlet prior the prior mean was set to be 0.0001, 0.001, 0.1, 1, 50, 100. Parameters α=c=1 were fixed. For each parameter was set four MCMC chains, and 2'500.000 generations and the analysis was run until the average standard deviation of split frequencies < 0.01. For the default exponential prior was used the same values of prior mean. Then, two nodes were extracted to see the behavour, one with a large branch and other one wihout large branchs.  To see the variation with respect to branch length was carried out the metric: Branch Score Difference (Kuhner & Felsenstein, 1994) taken exponential tree as reference topology and Robinson Foulds metric (Robinson & Foulds, 1981) to see if the topology changed. I simulated one tree with a large branch and another one without large branch using Seq-gen (Rambaut and Grassly, 1997).

Results and discussion


Branch lenght:
In general there is a pattern when Bt obtains values greater than 1, BSD decreases. In other words the difference in branch lengths between the consensus trees obtained using exponential and Dirichlet prior is becoming smaller when Bt increases (Fig 1).

Node 133 is similar to reference tree (exponential). Unlike 99 node which presents more variation respect the tree obtained under exponential prior, this can be explain by the large branch length (Fig 1).

Figure 1. Branch Score Difference (BSD) between each consensus tree of the exponential and Dirichlet prior for both nodes with differents values of Bt (prior mean).

For both nodes the exponential prior produces larger trees than compound Dirichlet priors for differents values of Bt. node 99 presents larger differences because exponential prior favors large branch (Zhang et al 2012) (Table 1).

Bt Expo Dir
Node 99 Node 133 Node 99 Node 133
1 0.3734763 0.05280543 0.134272 0.029316
0.1 0.386117 0.05436114 0.143912 0.029326
0.01 0.3879789 0.05413017 0.134548 0.029298
10 0.3100098 0.04519354 0.126194 0.028405
100 0.1683404 0.03025578 0.09396 0.022991
0.0001 0.3879235 0.05424801 0.145419 0.029751
50 0.2114827 0.03458418 0.105251 0.024711

Table 1. Tree length (TL) for each concensus topology under differents Bt values for exponential and Dirichlet priors.                                          

Regarding topology there isn't difference among trees obtained from exponential and Dirichlet prior (RF = 0).


Posterior distribution:
In the posterior exponential prior generated longer trees than Dirichlet prior (Fig 2).

Figure 2.  Effect of the exponential and Dirichlet prior on the posterior distribution of tree length.. Each colour represents a differente BT value: Red=1, pink= 0.1, blue claro= 0.01, green= 10, dark blue=100 , orange=0.0001 , gray=50.

In the complete topology, exponential prior shows the largest phylogenetic trees especially when Bt takes values less than 1, when Bt takes values greater than 1 decreases the size. The standard distribution is larger for Dir prior than exponential prior being greater for Bt values less than 1 (Fig 3).
Figure 3.  Standard Deviation (SD) for the tree total length from each Bt value.

In the total topology the behavour is the same like 99 and 133 node. BSD begins to decrease for Bt values greater than 10 but RF shows differences.
The values obtained are very close to the reference topology (tree simulated with and without large branch) using both priors and in Dir prior the variance was a little bit bigger than Exp prior, which was also found in 99 and 133 nodes.
Despite the difference between the results obtained from an empirical perspective and simulations, I think that this differences relate to the fact that in simulations all is perfect, all is in order and in biological systems is difficult to reach this point. I also found that there isn't influence in the posterior distribution from prior distribution in the Dirichlet prior and in exponential prior this does not happen. For this reasons I consider the Compound Dirichlet prior could be the solution to branch-length overestimation in bayesian inference because this produces more reasonable branch-length estimates than the default prior used in MrBayes and suggest to do more replicates and work with a known taxonomic group for which you know its diversification times to reach a best view. 


Brown JM, Hedtke SM, Lemmon AR, Lemmon EM. 2010. When trees grow too long: investigating the
causes of highly inaccurate Bayesian branch-length estimates. Syst Biol, 59:145-161.

Edgar, R.C. (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput
Nucleic Acids Res. 32(5):1792-1797

Kuhner, M. K., & Felsenstein, J. (1994). A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates. Molecular Biology and Evolution, 11(3), 459–468.

Marshall DC. 2010. Cryptic failure of partitioned Bayesian phylogenetic analyses: lost in the land of
long trees. Syst Biol, 59:108-117.

Rannala B, Zhu T, Yang Z. 2011. Tail paradox, partial identifiability and influential priors in Bayesian
branch length inference. Mol Biol Evol.

R Core Team. (2014). R: A Language and Environment for Statistical Computing. Vienna, Austria.

Robinson, D., & Foulds, L. R. (1981). Comparison of phylogenetic trees. Mathematical Biosciences, 53(1), 131–147.

Schliep, K. P. (2011). phangorn: Phylogenetic analysis in R. Bioinformatics, 27(4), 592–593.

Zhang C, Rannala B, Yang Z. 2012. Robustness of compound dirichlet prior for bayesian inference of branch lenghts. Syst Biol.

Summarizing branch lengths in phylogenies

Phylogenetic trees are a way to represent the evolutionary relationships of organisms; using them in biology is very broad. On numerous occasions phylogenetic analyzes  yield more than a phylogenetic reconstruction, the task of deciding how to represent all this variety of results can be complicated.

However, the overwhelming majority of attempts are carried out on trees branches whose length is despised. Take into account the length of the branches makes sense in statistical reconstruction methods, however the significance of these is difficult to determine. As an alternative to this deficiency one could take a measure of central tendency, the mean (average) or median of the branch lengths of all the alternative reconstructions. Felsenstein (2004) argues that one could choose as a parameter the median length of the branches on a tree consensus, however, there are problems when the clade of interest is not present in all the trees.

Take for example the case of support measures. Usually we see how the bootstrap values ​​and posterior probability are reported in the product trees ML and Bayesian Inference. But where exactly are these trees and branches lengths?

In Maximum Likelihood the result is a point estimate, the tree with the highest likelihood value. To take into account the confidence limits is generally used resampling matrix using bootstrap. The results are a measure of support from each of the internal nodes. Typically the Maximum Likelihood  tree  is reported and in it  is annotated the  node frequencies in the bootstrap replicates.

More complex Bayesian inference analysis in during each of the generations a different tree algorithm according to MCMC is sampled. What you get is a posterior distribution of trees rather than a single point estimate. This fact makes a representation to be obtained of the set of trees sampled.

The two software most used in  Bayesian inference in phylogeny are MRBAYES (Ronquist and Huelsenbeck,2003)  and Beast (Drummond et al, 2012) ; information of the run  can be summarized in a single tree. MRBAYES makes a majority rule consensus of  compatible trees  and if the ​​branch lengths values  have been saved, the average ​​branch lengths  and  their variance.. This tree can, in some cases, have not been found during the search conducted.

Beast, meanwhile, has a his own software to summarize the results of the trees: Treeannotator. The default option  for this utility is to produce a tree that maximizes the product of the posterior probability of bipartitions, called as maximum clade credibility tree. The branch length which  is taken may be the mean, median, common ancestror heigths  (hereafter ca) or maintain the target tree heigths. Helet and Buckaert (2013) show as there is no single methodology for summarizing  tree that maximizes the posterior probability of the clades and at the same time, has good branch lengths using different metrics.

Here, in a much simpler test than that used by them, I try to explore the differences between choosing different ways to summarize the trees of an  analysis of Bayesian inference., including topology (consensus, maximum clade credibility  ) and different branch length measurements  (mean, median, ca)


Two scenarios of evolution of sequences were explored; one under strict molecular clock and the other under uncorrelated relaxed molecular clock . For both scenaries , a tree were generated using RateEvolver (Ho et al, 2005), this software  produce a phylogeny of 9 terminals completely bifurcated , whose branch lengths measured as substitution  rates are modified according to the chosen molecular clock.
 From these trees,  10 replicates of 1000 nucleotides sequences  were generated in Seq-Gen (Rambaut and grass, 1997)  under the  JC model.
Phylogenetic reconstruction was carried out in Beast v1.8.3  following the molecular clock model under which the data were generated. 1500000 generations were run, 10% was used as a burn-in.
I summarized trees using two software:  Sumtrees (Sukumaran and Holder, 2015)  and Treeannotator  discarding the burn-in. Measurements of branch length   used were mean, median and ca. For topologies use MCCT and   majority  rule consensus.
It was calculated for each replicate the difference from the original tree, for this function treedist  of the phangorn (Schiliep, 2011) R package was used .

Results and Discussion

The results are shown in Figures 1 and 2.
Under a uncorrelated  relaxed molecular clock results are quite similar (Fig 1). The lowest values ​​of difference are using the majority rule consensus using either mean or median value as branch lengths. A similar value is found using  MCCT when ca is chosen to represent the branches.

Fig 1. Uncorrelated Relaxed  Clock. Results of mean Branch Score Difference for 10 replicates with respect to reference tree

The results obtained with strict molecular clock (Fig 2) results seem contradictory. The best distance values ​​(BSD) from the original tree are obtained using s maximum  calde credibility tree . A abrupt difference in the distance is observed compared to consensus trees. Given the kind of molecular clock it is expected that the results of consensus and MCCT would be very similar. One possible explanation could be a mistake in the software used, but  this seems unlikely given that  the results obtained under relaxed molecular clock seem fair.

Fig 2. Strict Clock. Results of mean Branch Score Difference for 10 replicates with respect to reference tree

With the results is not possible to draw a clear conclusion about what is the best way to summarize the phylogenetic reconstruction with branch length taken into account. The difference between choosing the mean or median does not seem to greatly affect the length of the branches. However, more work must be done due to ambiguous results reported here.


Drummond AJ, Suchard MA, Xie D & Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST 1.7 Molecular Biology And Evolution 29: 1969-197

Felsenstein, J., & Felenstein, J. (2004). Inferring phylogenies (Vol. 2). Sunderland: Sinauer Associates.

Heled, J., & Bouckaert, R. R. (2013). Looking for trees in the forest: summary tree from posterior samples. BMC evolutionary biology, 13(1), 1.

Ho, S. Y., Phillips, M. J., Drummond, A. J., & Cooper, A. (2005).
Accuracy of rate estimation using relaxed-clock models with a critical focus on the early metazoan radiation.
Molecular Biology and Evolution, 22(5), 1355-1363.

Rambaut, A., & Grass, N. C. (1997). Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees.Computer applications in the biosciences: CABIOS, 13(3), 235-238.

Ronquist, F. and J. P. Huelsenbeck. 2003. MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19:1572-1574. 

Schliep, K. P. (2011). phangorn: Phylogenetic analysis in R. Bioinformatics, 27(4), 592–593.

Sukumaran, J and MT Holder. SumTrees: Phylogenetic Tree Summarization. | | 4.0.0 (Jan 31 2015). Available at | | https://github.com/jeetsukumaran/DendroPy.

lunes, 28 de marzo de 2016

Optimization of continuous characters in parsimony

The node state in a phylogenetic tree can be easily calculate using discrete characters considering that the state is one option among a finite number of possibilities. However, calculate a node state using continuous characters is very complicated because the numbers of states are infinite (Goloboff et al, 2006). For this reason, continuous characters must be discretized. Farris (1970) proposes to treat continuous characters as additive characters. Several algorithms have been implemented (Goloboff, 1993 ; Thiele, 1993) to discretized continuous characters. Thus, I evaluated a Goloboff's algorithm and Thiele's algorithm by the Robinson-Foulds Metrics and Consistency Index (CI) expecting that the level of homoplasy decreases using continuous characters added to discrete characters. In order to, I simulated a tree of 25 tips (Figure 1) using geiger R library (Harmon et al, 2008), also simulated 3 different types of continuous characters and DNA sequence with 500 bp, model JC using seq-gen (Rambaut and Grassly, 2001). I recovery the tree in TNT (Goloboff, 2008) using only continuous characters, only discrete characters, continuous and discrete characters, different numbers (3, 10, 20) of states for thiele's algorithm and Goloboff algorithm, after I estimate the CI (Figure 2, Figure 3, Figure 4) and compare the trees with the original using Robinson-Foulds metric (Steel and Penny 1993) (Table 1).

Figure 1. Used simulated tree for generate continuous and discrete data.

Figure 2. Consistency Index for continuous (continuos) and discretes characters (discretos)

This could be because only three continuous characters were used. However, when more discrete continuous characters were used the degree of homoplasy decreased (Figure 3), this is consistent with the results obtained by Goloboff et al (2006).

Figure 3. Consistency Index for discrete characters (discretos) and continuous plus discrete characters (mix goloboff).

Figure 4. Consistency Index for total evidence using the Goloboff algorithm and thiele algorith with 3, 10 and 20 numbers of states.

Table 1. Robinson-Foulds Metric.

Goloboff algorithm yielded best results in the optimization way of characters. The problem with the Thiele’s algorithm is that always is subject to the subjectivity of the choice of the number of possible states. However, the Goloboff’s algorithm is not necessarily the best option for the reconstruction of phylogenies using total evidence. Robinson-Foulds metric (Table 1) showed that the smallest difference between the tree recovered and the original was obtained using the Thiele’s algorithm.

A copy of the scripts and used data can be found in https://github.com/dpabon/bio_comparada/tree/master/continuous_parsimony 


Farris, J., 1970. Methods for computing Wagner trees. Syst. Zool. 19,83–92. 
Goloboff, P., 1993b. Character optimization and calculation of tree lengths. Cladistics 9, 433–436.
Goloboff, Pablo A., Camilo I. Mattoni, and Andrés Sebastián Quinteros.' Continuous Characters Analyzed as Such’. Cladistics 22, no. 6 (December 2006): 589–601. doi:10.1111/j.1096-0031.2006.00122.x.
Goloboff, Pablo A., James S. Farris, and Kevin C. Nixon. ‘TNT, a Free Program for Phylogenetic Analysis’. Cladistics 24, no. 5 (1 October 2008): 774–86. doi:10.1111/j.1096-0031.2008.00217.x.

Harmon Luke J, Jason T Weir, Chad D Brock, Richard E Glor, and Wendell   Challenger. 2008. GEIGER: investigating evolutionary radiations. Bioinformatics 24:129-131.

Rambaut, A. and Grassly, N. C. (1997) Seq-Gen: An application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees. Comput. Appl. Biosci. 13: 235-238.

Steel M. A. and Penny P. (1993) Distributions of tree comparison metrics - some new results, Syst. Biol.,42(2), 126-141

Thiele, K., 1993. The Holy Grail of the perfect character: the cladistic
treatment of morphometric data. Cladistics 9, 275–304.

viernes, 18 de marzo de 2016

Statistical consistency of Maximum Likelihood and Maximum Parsimony

Statistical consistency is often one of the factors often cited to favor statistical estimation methods in phylogenetic reconstruction . It is defined as the porperty of estimator, for wich as data (characters)  tends to infinity, the inferred topology lead to  converge to the true: the real phylogeny (Truszkowski, & Goldman, 2015).

Debate about inconsistency  has a great story behind it. Since the proposal Felsenstein(1978) with a simple evolutionary model, much work has been done to evaluate the statistical properties of estimators in phylogeny under a set of different scenarios. Much of the effort has made to show that model based approaches like Maximum Likelihood (and Bayesian Inference) has advantages over Parsimony (Wheeler, 2011).

Here I test the statistical consistency of Maximum Likelihood  and two flavours of parsimony reconstruction in different sets of  topologies and models of nucleotide substitution,  with four taxa-trees  in a simulated experiment.

Materials and methods

To carry out the analysis of statistical consistency. First nucleotide sequences were simulated in the software Seq-Gen (Rambaut & Grass, 1997). in three different models of DNA evolution: GTR, and JC HKY. Also  each model was analyzed in four topologies (original topologies) whose length (q and p) differs as defined below:

  • Farris Zone.  q=0.6, p=0.1.Tips with long branches are sisters .
  • Felsenstein Zone. p=0.6, q= 0.1.  Tips with branches of different length are sisters.
  • Long Branches. p=q =0.6
  • Short Branches.p=q= 0.1
  • Central branch of topologies always take the value of q.

For maximum parsimony analysis two approaches were used: Static and dynamic homology (The static homology parsimony was conducted in TNT (Goloboff et. al, 2008). POY ( Varón et al, 2010) was the software used in the analysis of dynamic homology. Maximum Likelihood inferences were carried out in PhyML software (Guindon et al, 2010).

To asses differences between the inferred topologies and the originals, the symmetric difference (Steel & Penny, 1993) was used as implemented in RF.dist function ot the R package  phangorn (Schliep, 2011). The frequency of the true tree recovered  was the measure of consistency here used.

Results and Discussion.

Results show that  maximum parsimony analysis are consistent in three of the four areas evaluated or topologies. These are the Farris zone, short branches and  Long Branches. Maximum likelihood was equally consistent in three zones: Long branches, short branches and Felsenstein zone (Figs 1,2,3).
There is no dramatic difference between frequencies of recovered rigth trees by model.  

The tendency is to converge to the rigth tree in three of the four zones as the data avalaible increase, regardless the model of the data.

 Fig 1.  Frequency of recovered trees under the JC model. X axis show log of  base pairs. 
Fig.2.  Frequency of recovered trees under the HKY model. X axis show log of  base pairs.

Fig 3. Frequency of recovered trees under the GTR model. X axis show log of  base pairs. 

It can be seen as direct optimization leads to inconsistency much faster than does the traditional static parsimony under the JC model. These results confirm that, regardless of the concept of homology used,  parsimony  is inconsistent (Warnow, 2012) for the  Jukes Cantor  model  and more complex models in the Felsenstein zone.
 Statistical consistency is a desirable property in estimates of phylogenetic reconstruction (Wheeler, 2011). Under different configurations branches estimators behave similarly and all fail in at least some zones, even under simple models.
If the three (or two) approaches here assesed show a similar behavior (Failing to achieve perfect consistency in at least one zone), and perfect consistency is not reached one must choose an approach between them with freedom knowing that in some configuration of branch length it can be yield the rigth tree when the avalaible data is enougth.


Felsenstein, J. (1978). Cases in which parsimony or compatibility methods will be positively misleading. Systematic Biology, 27(4), 401–410.

Goloboff, P. A., Farris, J. S., & Nixon, K. C. (2008). TNT, a free program for phylogenetic analysis.Cladistics, 24(5), 774-786.

Guindon, S., Dufayard, J. F., Lefort, V., Anisimova, M., Hordijk, W., & Gascuel, O. (2010). New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Systematic biology, 59(3), 307-321

Rambaut, A., & Grass, N. C. (1997). Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees.Computer applications in the biosciences: CABIOS, 13(3), 235-238.

Schliep, K. P. (2011). phangorn: Phylogenetic analysis in R. Bioinformatics, 27(4), 592–593.

Steel, M. A., Hendy, M. D., & Penny, D. (1993). Parsimony can be consistent! Systematic Biology, 581–587.

Truszkowski, J., & Goldman, N. (2015). Maximum likelihood phylogenetic inference is consistent on multiple sequence alignments, with or without gaps.Systematic biology, syv089.

Varón, A., L. S. Vinh, W. C. Wheeler. 2010. POY version 4: phylogenetic analysis using dynamic homologies. Cladistics, 26:72-85.

Warnow T. Standard maximum likelihood analyses of alignments with gaps can be statistically inconsistent. PLOS Currents Tree of Life. 2012 Mar 12 . Edition 1. doi: 10.1371/currents.RRN1308.

Wheeler, W. C. Comparison of Optimality Criteria. Systematics: A Course of Lectures, 2011.269-287.

ML and IB consistency

Statistical consistency is the capacity of the estimator that be true when the lenght of data tend to infinite. Felsenstein (1978) show that phylogenetic reconstruction using parsimony is incosistent when the branch lengths have certain proportions (Figure 1) this proportions have the name of Felsentein zone. However Farris (1999) show that phylogenetic reconstruction using Maximum Likelihood (ML) is inconsistent also when branch lengths have certain proportions (Figure 2) this type of zone have the name of Farris zone.
Contrary to the previous phylogenetics reconstructions techniques, bayesian analysis (IB) is consistency without import the zone or branch lengths relations (Steel, 2013).

Figure 1. Felsenstein zone

Figure 2. Farris zone

However i checked experimentally if is true the consistency of IB and ML. For this i use four type of unrooted trees with 4 tips (farris zone, felsenstein zone, "long branches", "small branches"). I choose three sizes base pairs (100, 1000, 10000) and for each size i generated 25 sequences of DNA based in the tree and model JC . I recover the trees using Mrbayes v 3.2.6 compiled to MPI and Phyml 20160303 also compiled to MPI, and i using the JC model. I compared each tree vs original tree using Robinson Fould metric (Steel & Penny, 1993). I define the consistency as: Number of trees that Robinson Fould metric = 0/Number of total trees for each type of tree and size of bases pairs.  
In the Farris zone (Figure 3) i need sample more sizes of DNA because only three sizes show that ML is consistent in the Farris zone and this is false, IB was consistent in this zone. In the Felsentein zone (Figure 4), ML and IB were consistent. In "long branches"(Figure 5) and "small branches" (Figure 6) ML and IB were consistent. For this is true that IB is consistency independently of the zone or branches lenght proportions. 

Figure 3. Statistical consistency in Farris zone. ML = Maximum Likelihood, IB =Bayesian Inference

Figure 4. Statistical consistency in Felsenstein zone. ML = Maximum Likelihood, IB =Bayesian Inference

Figure 5. Statistical consistency for "Long branches". ML = Maximum Likelihood, IB =Bayesian Inference

I programmed some useful functions in R for the threes.
If you want run all analysis in your computer you can download a copy of the repository and follow the instructions.


Farris, J. (1999). Likelihood and Inconsistency. Cladistics 15, 199–204

Felsenstein, J. (1978). Cases in which parsimony or compatibility methods will be positively misleading. Systematic Biology, 27(4), 401–410.

Steel M. A. and Penny P. (1993) Distributions of tree comparison metrics - some new results, Syst. Biol.,42(2), 126-141

Steel, M. (2013). Consistency of Bayesian inference of resolved phylogenetic trees. Journal of theoretical biology 336: 246–49.