Lab 2: Bayesian Phylogenetics with MrBayes

Naxos Summer School 2026 · Phylogenetic Methods in Historical Linguistics

Nothing to install

This page runs R inside your browser. No R to install, no MrBayes, no account, and nothing you write is sent anywhere. (R itself and its packages are fetched from a public CDN on first load.) Edit any cell and press Run Code.

The page takes about a minute to set itself up. It gives no sign of progress while it does so, and then everything appears at once. This is normal.

The analysis has already been run

A Bayesian analysis does not finish inside a lecture slot: the one below ran for 6,453,000 generations, about twenty minutes. MrBayes has therefore been run in advance, and what we have is its output. That is what one would spend one’s time on in any case. Nobody’s professional life consists of watching a chain run. It consists of deciding whether to believe what came out of it.

That is what this lab is about. We are handed the output of a Bayesian phylogenetic analysis, and we have to decide whether to trust it.

What went in

The binary cognate matrix from Lab 1: 13 Indo-European languages by 1,030 variable characters. Each character is one cognate set, coded 1 if the language has a reflex of it and 0 if not.

lset  coding=noabsencesites;                  [ascertainment correction, cf. Lecture 4]
lset  rates=gamma ngammacat=4;                [Gamma rate variation]

prset brlenspr=clock:birthdeath;              [a ROOTED tree: birth-death process prior]
prset speciationpr=exponential(10);
prset extinctionpr=beta(1,1);
prset sampleprob=0.03;                        [13 of ~450 IE varieties]
prset clockvarpr=igr;                         [RELAXED clock: rates vary across lineages]
prset igrvarpr=exponential(10);

calibrate Latin         = fixed(2100);        [<-- the ancient languages ARE the calibration]
calibrate AncientGreek  = fixed(2450);
calibrate VedicSanskrit = fixed(3200);
prset nodeagepr=calibrated;
prset treeagepr=offsetexponential(3200,6000);
prset clockratepr=lognormal(-8.5,1.5);

mcmcp ngen=10000000 samplefreq=2500 nruns=2 nchains=4;
mcmc;                                         [stopped by hand at 6,453,000]
sumt relburnin=no burnin=645;                 [2,582 samples per run, 25% discarded]
sump relburnin=no burnin=645;

Three things in that block do the real work.

brlenspr=clock:birthdeath makes the tree rooted. Lecture 4 showed that a reversible substitution model cannot find the root, since the likelihood is identical wherever we put it. Information from outside the substitution model is required, and a clock is one way to supply it. The birth-death process is a generative model of how lineages split and die out, and it gives us a prior over rooted, time-scaled trees.

The three calibrations turn substitutions into years. Latin, Ancient Greek and Vedic Sanskrit are not contemporaries of the modern languages, and a clock model would otherwise force them to the present, which is simply false. Telling MrBayes how old they are corrects this, and as a bonus it anchors the whole tree in absolute time. The ancient languages are thus not merely data. They are the clock.

clockvarpr=igr relaxes the clock. A strict clock would insist that every lineage changes at the same rate. It does not: Icelandic is famously conservative, English is not.

Two settings are constantly got wrong:

  • nruns=2 gives two independent runs, and convergence is judged on these. nchains=4 gives four Metropolis-coupled (“heated”) chains within each run. They help the cold chain escape local optima, they are not extra samples, and they never produce a second .p file.
  • relburnin=no is not optional. MrBayes defaults to relburnin=yes with burninfrac=0.25, so a bare burnin=1000 is silently ignored.

1. Did it converge?

This step should not be skipped. A posterior from a chain that has not converged is not a posterior. It is a pile of numbers with a mean.

There is no need to compute PSRF ourselves. MrBayes has already done so and written it to .pstat.

TH is tree height and TL tree length, both in expected substitutions. clockrate converts substitutions to years. net_speciation and relative_extinction are the birth-death parameters. alpha is the Gamma shape for rate variation, and igrvar the variance of the relaxed clock.

PSRF sits at 1.00 across the board (the largest is 1.009), and every ESS clears the 200 threshold, though the smallest, 246 on the state frequencies, does not clear it by much. MrBayes also reports an average standard deviation of split frequencies of 0.003, well under the 0.01 rule of thumb. Three independent signals thus agree that this run is sound.

Clock parameters mix slowly, and PSRF can look perfect while the effective sample size is tiny: an earlier attempt at 2,000,000 generations gave PSRF = 1.000 and ASDSF = 0.006, which looks fine, alongside an ESS of only 56 on tree height and clock rate. That is why this run was carried to 6,453,000, and why the ESS column matters as much as the PSRF one. Both need to be checked.

2. Look at the chains anyway

The first quarter of each chain is the burn-in and has been discarded, so both panels show only what we would actually keep. What we want to see is a fuzzy caterpillar with no trend across it, two runs overlapping so completely that they cannot be told apart, and two densities lying on top of each other.

Failure looks quite different. Each chain looks perfectly smooth on its own, but the two densities sit apart as separate humps. This is the trap: a single chain can look immaculate and still be wrong, and the problem cannot be detected with one run. That is what nruns=2 is for.

3. The tree, rooted and in years

The tools here are the real ones

Everything below is stock ggtree and treeio, i.e. the packages this field actually uses. They are Bioconductor, and webR’s repository is CRAN-only, but both are pure R, so we simply unpacked an installed build (see the setup cell). Any of this code can be copied into RStudio at home, where it runs unchanged.

treeio::read.mrbayes() matters more than it looks. MrBayes stores clade support in the consensus tree as a NEXUS comment, [&prob=0.73,…], and plain ape::read.nexus() discards comments, so the support silently vanishes. read.mrbayes() keeps it, and returns a treedata object carrying both the tree and its annotations.

Read the support numbers before reading the topology. The five branches, Germanic, Romance, Slavic, Hellenic and Indic, all come out at posterior 1.00, which is why they carry no label. Every number that is printed sits on a node in the deep part of the tree, and none of them reaches 0.75.

The chain is thus entirely sure which languages group into which branch, and close to clueless about how those branches relate to one another. That contrast is not a defect of the analysis. It is the honest shape of the evidence, and it is why the Indo-European debate is still a debate.

But a consensus tree is a lie of omission

That single tree is a summary. It shows one topology and one set of dates, and it draws them with equally confident lines whether the posterior was certain or hopelessly torn.

Let us therefore look at the posterior itself. A DensiTree overlays every sampled tree, one faintly on top of another. Where the trees agree, the lines pile up and go dark. Where they disagree, the ink spreads out into a haze.

The uncertainty can now be seen, rather than merely asserted.

Look at the tips: sharp, single lines. We know where the languages are, and the three ancient ones stop dead at their calibrated dates, exactly as instructed.

Look at the shallow nodes, i.e. Germanic, Romance and Slavic: tight, dark, decisive. The posterior is practically unanimous.

Now look at the deep nodes, near the root. The ink is smeared across thousands of years. That haze is not a rendering artefact. It is the answer, and it is the same uncertainty that the consensus tree drew as a single confident black line.

A consensus tree tells us what the posterior thinks. A DensiTree tells us how strongly it thinks it.

The numbers on that axis now mean something. The ancient languages sit back in time where they belong, the modern ones at zero, and every internal node has an age.

All five branches are recovered. Now let us argue with the dates. Slavic breaking up only about 700 years ago is very young, and Germanic at about 1,360 BP is plausible for the West/North split but not for Proto-Germanic itself. Are those failures of the model, or artefacts of having only two or three languages per branch? What would we need to check?

4. The number everyone actually wants

The root of this tree is Proto-Indo-European, and the whole of Lecture 5’s debate is about how old it is. We have a posterior, so let us look at all of it, and not only at the point estimate.

Read that plot before reading the next sentence

Our median is about 6,000 BP, which is the steppe answer. Anyone stopping here would write a paper saying so, and would be able to cite a Bayesian analysis with converged chains and excellent diagnostics in support of it.

Now look at the 95% credible interval, which runs from roughly 4,700 to 8,200 BP. It still contains Heggarty’s 8,100, though only just, and it falls short of Gray and Atkinson’s 8,700.

The point estimate thus appears to settle the debate, while the interval settles nothing. Thirteen languages and 170 concepts cannot distinguish the steppe hypothesis from the Anatolian one, and the honest report of this analysis is that it does not.

This is the whole course in one figure. A result need not be believed merely because it has a number attached, and a median may not be quoted without its interval.

5. And what if we had not used a clock?

Without clock:birthdeath we get an unrooted, undated tree, and brlenspr reverts to a prior on branch lengths directly. The default there is gammadir(1,0.1,1,1), the compound Dirichlet–Gamma prior.

On this data, that default does not converge. Not slowly, but not at all, even after two million generations:

The two runs settle in different places and stay there. Tree length differs about fivefold, α by more than an order of magnitude, and PSRF sits above 2 and 3. The recommended, modern, default choice is thus the one that failed.

The prior is therefore not a technicality tucked away in the methods section. It is the difference between an answer and a mirage.

Exercises

1. Would more data narrow that interval?

The credible interval on PIE spans about 3,500 years. Where would you spend effort to narrow it: on more languages, more concepts, or more calibrations? Consider what each of them actually constrains. Which of the three directly informs clockrate, the parameter that converts substitutions into years?

2. Break a calibration

In data/precomputed/iecor_bdp.nex, Vedic Sanskrit is fixed at 3,200 BP. Suppose the true figure were 3,800. Without re-running anything, reason it through: does the PIE estimate get older or younger, and roughly by how much? Then argue why a fixed calibration is dishonest, and what you would use instead.

3. The clock you did not question

We used a relaxed clock (igr). Try to justify a strict one for Indo-European, then try to demolish your own justification. Icelandic and English have been separated for about 1,100 years; do they look equally changed? Which assumption is doing more work in the date above, the clock model or the calibrations?

4. Break the burn-in

In §2, change burnin <- round(n * 0.25) to round(n * 0.9), and then to 0. Re-run §2 and §4. What happens to the credible interval on PIE, and which of the two mistakes would do more damage in a published paper?