The Subulate Palpomere

Many years ago I dressed up for Halloween with a piece of paper taped to my rear with “Pr” on it, in one hand holding a small crowbar, and in the other an engagement ring in a small box (see below). The “Pr” was a posterior probability, the crowbar a prior (a “pryer”), and the engagement ring a proposal mechanism – that was, I was dressed up as a Bayesian analysis.

You might think that would suggest I love Bayesian analyses, but whenever I think of setting up a Bayesian analysis to infer aspects of phylogeny, my heart sinks. I hate wandering into the gauntlet of choices one has to make – it’s as bad as filling out tax forms, or perhaps even worse, as for tax forms there is sufficient documentation available to help me, and I have pieces of paper that record what values I should use.

Among the model-based methods for inferring phylogenetic trees, the two most popular are Maximum Likelihood Estimation (MLE) and Bayesian approaches. The former use what is sometimes called a “frequentist” approach to probabilities, and the latter a Bayesian approach. Bayesian analyses tend to include a richer and more complete model of the evolutionary process, and they require a specification of our prior beliefs about it. I am someone used to the simplicity of MLE approaches, and the relative lack of choices to be made. Diving into Bayesian analysis requires me to face my ignorance about the evolutionary processes that generated the beetles’ diversity, and the processes we have used to sample the beetles.

I am currently attempting to conduct a Bayesian analysis of some genomic data in order to estimate when some lineages of beetles diverged one from another. This is my first attempt to do such an analysis to date phylogenetic splits, and my goodness I am finding it challenging. I decided to expose my confusions and queries to the world, and hope that a kind soul who knows much more about Bayesian analysis than me answers me call. I’ve already received a bit of feedback from Sean Harrington to an earlier set of questions (thank you, Sean!). I hope that in addition to answering my questions, this current effort might help someone else in the same predicament. My goal is to update this post with notes as I get feedback about how best to make the choices needed.

For me the most difficult part of setting up a Bayesian analysis is specifying the priors used. A prior is a statement of the relative probabilities of the potential values of various components of the model, with those probabilities based upon our prior knowledge (rather than the data at hand). However, so little is known about the beetles and the evolutionary process that to specify these priors seems foolhardy; at times I feel as if I need a magic 8-ball or tea leaves or access to those rarefied few who know the secret chants. But Bayesian analyses require such choices to be made, and so I have to make them, and I would rather not trust magic 8-balls or tasseography.

Bayesian analyses can take a very long time (especially if the model used is complex), and that in itself leads to a desire to make good decisions about these choices from the start, as it is time-costly to start an analysis, only to find out several weeks later that mistakes were made.

There are several programs available for conducting Bayesian analyses of phylogenies, including the venerable MrBayes and its descendant RevBayes, BEAST2, and PhyloBayes. My initial explorations suggest that RevBayes and PhyloBayes might require more computer processing power than I have easily available, and so I have decided to use the popular BEAST2.

For a long time there was very little guidance to the user regarding how the many choices involved in setting up a Bayesian analysis should be made, except for the occasional post here or there on a discussion forum. There is one extremely useful document by Tracy Heath, called Divergence Time Estimation using BEAST v2.∗ Dating Species Divergences with the Fossilized Birth-Death Process. However, the ever-evolving suite of models used in the BEAST2 often requires a new set of choices when a new version comes out, and that document is now out of date. More recently, the BEAST2 developers have added an extremely useful “Help Me Choose” site. However, many of the choices I need to make aren’t discussed there yet, especially for some of the newer additions to BEAST2.

Before I go through what choices I made (and where my confusions are), I’ll give some background about the data set I happen to be dealing with at the moment.

The taxonomic group and the data

The data matrix consists of about 500,000 amino acids from many genes for 46 species of beetles (see gratuitous beetle picture below). Of those 46 species, four are outgroups (that is, they do not belong to the group I am focusing on, and are included to help root my study group). The age of the clade represented by the entire 46 species is at least 99 million years, based upon some undescribed fossils; the full clade has about 5500 described species, and likely at least that many again undescribed ones. The ingroup (the primary study group) of 42 species has about 1300 described species, but there are likely at least 2,000 species currently living. In addition, there are three well-documented fossils described in the ingroup, all of which are between 34-48 million years old (the fossils haven’t been more accurately dated than that). There are no estimates about the age of the ingroup beyond the realization that the group is at least 34 million years old. I did not partition the 500,000 sites as initial tests suggested the analysis would take much too long if I did.

Setting up an analysis for BEAST2 in BEAUTi

The following analysis was set up for BEAST2 version 2.7.3. In addition to the core BEAST2 packages, I also installed bModelTest and OBAMA. Those packages allowed me to do Bayesian amino-acid model averaging. Here’s what I did:

• Open BEAUTi
• Choose File > Manage Packages
• In the window that appears, select bModelTest, and press Install/Upgrade
• Select OBAMA, and press Install/Upgrade
• Quit BEAUTi, and restart it

I prepared my NEXUS file by adding three new taxa in Mesquite, one for each of the fossils. Thus, in total the matrix has 49 taxa in it (four outgroups, 42 ingroups, and three ingroup fossils containing no DNA sequence data).

I then created three taxon sets in Mesquite; each taxon set included one of the fossils, and the remaining terminal taxa that form the smallest clade to which the fossil definitively belongs. We know that each of those three clades is at least 34 million years old. After saving the file, I opened up the file in a text editor, stripped out extra NEXUS file commands, and replaced the numbers in the TAXSET commands with the full taxon names. (I’ve written a new Mesquite module to Export files for BEAUTi to take care of these things automatically, and that will come out in the next release of Mesquite.)

I then loaded my NEXUS file containing the data into BEAUTi by choosing File > Import Alignment.

Here are the options I chose in BEAUTi, panel by panel. My questions/confusions are shown in blue.

Update: the team at the Centre of Computational Evolution at University of Auckland, keepers of the BEAST, were kind enough to prepare a response. The response was a joint effort by Kylie Chen, Alexei Drummond, Remco Bouckaert, and Walter Xie. Their full response is in their comment, below, but I have also inserted responses to individual questions within the blog post, for easy reading. Where I felt a response was appropriate, I have added them in green. These include addition information about what I chose in response to their answers.

---

Tip Dates Panel

I feel pretty comfortable about my choices for this panel. In particular, I

• Checked “Use Tip Dates”
• Switched to Dates specified numerically as year Before the present
• For the three fossils (Bembidion_christelae, Bembidion_alekseevi, and Bembidion_bukeisi), changed their “Data (raw value)” to 41. Each of these fossils is from Baltic amber, which is 34–48 million years before the present. I choose 41 as that is the center of that range.

Here’s what the the lower part of the Tip Dates panel then looked like:

(1) Bayesian gurus: are these good choices for the Tip Dates panel?

Response: Yes. Please note that if you sample uncertainty in the age of the fossil taxa then the 41 million you specify here is just a starting value and nothing more.

---

Site Model Panel

In this panel I:

• Chose OBAMA Bayesian Aminoacid Model Averaging
• Deselected the options that appeared to be specific to other types of sequences, and not to metazoan nuclear protein coding genes

This is what the Site Model Panel then looked like:

(2) Is it OK to leave Mutation Rate as it is?

Response: Yes. There is only one partition and therefore the partition-specific mutation rate should not be estimated or it will be non-identifiable with the clock rate.

(3) Are the seven models selected reasonable ones to include for an analysis of nuclear protein-coding genes in an animal?

Response: Yes, these are reasonable selections. The unselected models are specific for virus/mitochondria/reverse transcriptase. We note that keeping all models selected would also be acceptable as the data should converge on the better fitting model.

---

Clock Model Panel

In the Clock Model Panel I selected Optimized Relaxed Clock, and didn’t change anything else, so this is what the panel looked like:

(4) Should I leave Clock rate as it is?

Response: Yes. Estimating the clock rate is appropriate because you have calibration information in the form of fossil taxa.

---

Priors Panel

And now, to the tough part, the priors. Here’s what the panel looked like at the start:

For the Tree.t prior, I choose Fossilized Birth Death Model. I then opened up that prior, and changed only one value, that for Rho. According to Tracy’s tutorial, rho is the probability of sampling a tip in the present. Presuming that there are about 10,000 living species in this whole group, and that we have sampled 46 of them, then rho should be 0.0046. I used the value 0.005. I set Origin to 120, as that is older than the oldest fossil. The Tree.t prior information thus looked like this:

(5) Is this reasonable? Should I change anything else, especially in the face of my almost complete ignorance of the evolutionary process in this group?

Response: The numbers entered as initial values do not affect the meaning of the model, however tweaking them may be needed to ensure the initial state is valid and can improve the time taken to achieve convergence (get through burn-in) during MCMC. You should choose a prior for your parameters centred around your known value. For example, a (Beta) prior on sampling proportion centred around 5E-4. The starting value can then be any value chosen from the prior distribution.

My response: Because I am now sampling only the ingroup (see the next question), I’ve reduced the initial values to reflect this smaller group, with 42 species of a possible 2,000 species sampled, and a younger age for the origin. Here are my revised choices:

(6) One thing I really struggle with here. Is it reasonable to use an FBD model when we know the sampling within the whole clade is not uniform? I sampled the ingroup much more heavily than the outgroup. I sampled only 4 species in the outgroups (of over 4,200 known species), but I sampled 42 ingroup species (of 1,300 known species). How can that differential sampling intensity be considered?

Response: The currently available FBD model does not account for variation in sampling rates among lineages. To minimise the effect of variation in sampling intensity we suggest two strategies: (a) remove the outgroup and run analyses using only the ingroup, (b) sample the outgroup as heavily as the ingroup.
In future we anticipate that extending the multistate birth death package MSBD to handle fossil data would provide a good solution to this problem  https://taming-the-beast.org/tutorials/MSBD-tutorial/.

My response: I can’t sample the outgroup as heavily as the ingroup (that would require a huge amount of time and money to do the genomic sequencing), which means I will choose (a), removing the outgroup. That’s OK as I don’t really care about the dating of the outgroup, and all of my fossils are within the ingroup anyway. I also have a good idea as to where the ingroup is rooted based upon other analyses. Until the models can cope with differential sampling, this seems like the best option.

I wonder whether or not the FBD model is really the one I should be using anyway. The sampling of the ingroup was not done randomly; I very specifically chose to sample one or two species from each of the major lineages, so it was much more dispersed than equiprobable sampling.

Update February 2024: Alas, as it turns out running the FBD model in BEAST was not a good idea. For details, and an update on the story of this analysis, see my more recent post, Maddison’s Law of Phylogenetic Analysis.

As I could find no guidance and no reason to change the OBAMA priors, I left them the way they were:

(7) Are these reasonable OBAMA priors?

Response: Yes, these defaults are justified in the OBAMA paper:

• OBAMA proportion invariant: Beta(1,4) has a mean of 0.2, and favours lower proportion invariable sites, but still allows large proportions.
• OBAMA gamma shape: any shape value <0.1 leads to one or more categories effectively being zero, which is what the proportion invariable category already models if your data contains a significant number of invariable sites. The OBAMA paper has more details/graphs on why the 0.1 cutoff is a practical choice.

Similarly, I could find no guidance and no reason to change the Optimized Relax Clock priors, and so I left them the way they were:

(8) Are these reasonable ORC priors?

Response:

• ORCRates: Yes. This is the prior distribution of the branch rates under the relaxed clock model, which are assumed to be drawn from a lognormal distribution with a mean of 1 (in real space) and a standard deviation of ORCsigma, below.
• ORCsigma: Yes. This default prior is reasonable for the general case. A standard deviation of 0.1 – 0.6 allows the branch rates to vary somewhat, but if this term exceeds 0.8, then this indicates the data are non-clock like in which case many vastly different trees may explain the data. This prior is centred around a realistic range of values.
• ORCucldMean: No, this is not a good prior choice. This is the clock rate prior. It should be informed by known related systems, and in the case of this beetle dataset it will be expressed in units of substitutions per site per million years. Since you have a calibrated analysis, you could assume a relatively uninformative broad log-normal prior on this parameter quite safely.

For the diversificationRateFBD.t prior, Tracy’s says in her FBD tutorial “Generally, we think that this value is fairly small, particularly since we have few extant species and many fossils. Therefore, an exponential distribution is a reasonable prior for this parameter as it places the highest probability on zero”. In my case there are lots of extant species, and few fossils. I have no idea what to choose, but I was advised “even in large clades, this is probably a relatively small value and an exponential distribution should be fine. If the actual rate is higher, an exponential doesn’t truncate higher values, just puts less weight on them.” I thus chose the default exponential:

(9) For a group like this with many species and few fossils, is a default exponential a reasonable choice for this prior?

Response: To assess whether the choice of prior is appropriate, you can look at how the posterior changes compared to the prior, or conduct a sensitivity analysis using different priors for the diversification rate.

For the originFBD prior, I was advised that a lognormal with a wide variance with peak where I guess it should be would be reasonable. So I chose that, setting the lower bound to 99 (as that is the age of the oldest fossils within the entire group), the upper bound to 250 as, based on the dating that has been on on Coleoptera diversification, 250 is surely beyond the maximum age of the entire clade (outgroup+ingroup); I set the initial age to 110.

I then set the offset to 99, the M to 120, and the S to 1. That gave a curve that feels reasonable for the origin time of the whole clade. That’s not based on clear evidence, but hopefully the curve is flat enough to be acceptable.

(10) Should I really be using a lognormal here? A uniform would also be possible, set between 99 and 250, but it does seem better to put more of the prior distribution at lower values. But how do I choose details of the lognormal? Does it matter?

Response: Hard boundaries reflect 100% certainty the value will not exceed those boundaries. We suggest using a smoother left skewed distribution with (99, 250) in the 95% interval rather than hard boundaries.

My response: I modified the log normal so that it doesn’t have hard boundaries. Also, because I am now including only the ingroup, the values have changed: the peak of the distribution is now much lower. Here’s what I now have:

I’m rather lost about what to choose for the samplingProportionFBD. Tracy’s tutorial says “The sampling proportion is the probability of observing a lineage as a fossil before that lineage goes extinct.” If there are about 10,000 living species. We have only three fossils. So, as the sampling proportion is very small, I chose an exponential distribution, with an initial value of 0.0005, and a mean of 0.01:

(11) Well? Reasonable choices?

Response: Sampling proportion FBD = the probability of sampling prior to death = sampling rate/(sampling rate + death rate). An uninformative prior on this parameter would be a uniform prior between 0 and 1. An informative prior would be some form of Beta prior that was biassed towards the 0 end, since you only have 3 fossils in your analysis compared to many more unsampled fossil species.

My response: The sampling proportion for the ingroup is 3 fossils in perhaps 2,000 species, so about 0.0015; I used that as my initial value and the center for the prior distribution. In the response to question 5, the Bayesian gurus suggested “a (Beta) prior on sampling proportion centred around 5E-4”. Here’s what I now have (but with a different center, around 0.0015, because of the removal of outgroups):

For turnoverFBD, I used a uniform 0 to 1 distribution, as in Tracy’s tutorial:

Each of the three taxon sets (one for each fossil) contain the set of species that form the smallest clade that we confidently believe contains that fossil (the fossil is also include in the set). Following Tracy’s tutorial, I assigned the prior for each of these taxon sets to follow a uniform distribution with the age range of the fossil contained in that taxon set. So all three taxon set priors then looked like this:

I think we are done with the priors! Here’s an overview of the BEAUTi Prior panel at the end:

At this point I feel about the same as I feel when I have got to the last page of my tax forms…

The MCMC Panel is more in my comfort zone, and I changed only two things: (1) chain length (which I increased to 100M, just in case), and (2) number of initialization attempts (to 100).

---

So I saved the file in BEAUTi, and tried to execute it in BEAST 2.7.3 with the following command:

[path to beast] -seed 22 -working -threads 8 -instances 8 [path to xml file]

Everything started up fine, but then I get a notice about a failure to find an initial starting point:

===============================================================================
Start likelihood: -Infinity after 100 initialisation attempts
P(posterior) = -Infinity (was -Infinity)
  P(prior) = -Infinity (was -Infinity)
    P(FBD.t:46Taxa_Occ66_AA_AllLoci) = -297.30215228670386 (was -297.30215228670386)
    P(OBAMA_PropInvariablePrior.s:46Taxa_Occ66_AA_AllLoci) = 1.0702128141464131 (was 1.0702128141464131)
    P(OBAMA_freqsPrior.s:46Taxa_Occ66_AA_AllLoci) = 53.71197185321927 (was 53.71197185321927)
    P(OBAMA_GammaShapePrior.s:46Taxa_Occ66_AA_AllLoci) = -1.0 (was -1.0)
    P(ORCRatePriorDistribution.c:46Taxa_Occ66_AA_AllLoci) = -410.92248229518384 (was -410.92248229518384)
    P(ORCsigmaPrior.c:46Taxa_Occ66_AA_AllLoci) = 1.3628558876856076 (was 1.3628558876856076)
    P(ORCucldMeanPrior.c:46Taxa_Occ66_AA_AllLoci) = -2.7232296703330143 (was -2.7232296703330143)
    P(diversificationRatePriorFBD.t:46Taxa_Occ66_AA_AllLoci) = -1.0 (was -1.0)
    P(originPriorFBD.t:46Taxa_Occ66_AA_AllLoci) = -5.102121215690426 (was -5.102121215690426)
    P(samplingProportionPriorFBD.t:46Taxa_Occ66_AA_AllLoci) = 4.555170185988091 (was 4.555170185988091)
    P(turnoverPriorFBD.t:46Taxa_Occ66_AA_AllLoci) = 0.0 (was 0.0)
    P(Eupetedromus.prior) = -Infinity (was -Infinity)
    P(Ocydromus_SuperSeries.prior) = NaN (was NaN)  **
    P(Philochthus.prior) = NaN (was NaN)  **
  P(likelihood) = NaN (was NaN)  **
    P(treeLikelihood.46Taxa_Occ66_AA_AllLoci) = NaN (was NaN)  **

java.lang.RuntimeException: Could not find a proper state to initialise. Perhaps try another seed.
See http://www.beast2.org/2018/07/04/fatal-errors.html for other possible solutions.
  at beast.base.inference.MCMC.run(Unknown Source)
  at beastfx.app.beast.BeastMCMC.run(Unknown Source)
  at beastfx.app.beast.BeastMain.main(Unknown Source)
  at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
  at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)
  at java.base/jdk.internal.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)
  at java.base/java.lang.reflect.Method.invoke(Unknown Source)
  at beast.pkgmgmt.launcher.BeastLauncher.run(Unknown Source)
  at beast.pkgmgmt.launcher.BeastLauncher.main(Unknown Source)

ARGGGHHHH!

I presume this is because there is some contradiction in the priors I chose. I’ll work on that issue (but if anyone has suggestions, they would be welcome!).

Response: Finally, the error message comes from your set up of the fossil calibrations. There are two logically distinct things that you want to specify for each fossil:

• Age uncertainty: specify the geological age uncertainty of the fossil.
• Topological constraint: specify the extant taxa that this fossil should be grouped with.

These two things need to be done with two different priors in BEAUti. So for three fossils you will specify 6 priors. For each of the fossils you will need specify two things:

• An ‘Sampled ancestor MRCA Prior’ containing only the fossil taxon (with tipsonly checked) and with a uniform distribution describing the geological age range of the fossil. This prior option is available in BEAUti with the latest version of the sampled-ancestors package when adding priors in the priors panel.
• An ‘MRCA Prior’ containing the fossil and the related extant taxa, with tipsonly unchecked and monophyletic checked. There is no age distribution associated with this prior, it is just used to maintain the topological constraint that the fossil must stay within the extant group it is associated with.

My response: At first I was a bit confused by this. It appears that the MRCA Prior is the one that I get by default in the BEAUTi interface when I defined the taxon sets in the NEXUS file that I loaded into BEAUTi to begin with. That is, the three priors that appear in BEAUTi (see below) based upon my pre-defined taxon sets are MRCA Priors; I can tell this not because it is indicated in the BEAUTi interface, but because that is what they are called in the XML file that BEAUTi saves.

To add the Sampled ancestor MRCA Prior, touched on the +Add Priors button at the bottom of the Priors panel, and choose Sampled Ancestors MRCA Prior from the dialog box that appears. (I get this only with the latest version of the SA package.) I then defined the taxon set to include only the fossil for that group, and set the prior to have a uniform distribution with the lower and upper bounds indicating the uncertainty in the fossil age.

For example, for one of the fossils, the following two priors were present, one for delimiting age uncertainty, and one for defining the clade:

The first one is the Sampled Ancestors MRCA Prior; the second one is the MRCA Prior.

There’s one more change I had to make. Because I removed all of the outgroups, I needed to enforce the rooting of the ingroup. Fortunately, from likelihood analyses, I am pretty confident that I know where the root of the ingroup goes: there is one species included among the 42 that is the sister group to all the rest. So I made a taxon set that includes all of the other 41 taxa, and made an MRCA Prior for it:

I made these changes, and ran it again. And now it works!

Any additional suggestions for improving the choices I made in priors and elsewhere would be most appreciated. You can leave a comment here, or you can email me at david.maddison@science.oregonstate.edu. Thanks!

Thanks so much, Kylie Chen, Alexei Drummond, Remco Bouckaert, and Walter Xie, for responding with answers to my queries!

Update February 2024: This analysis proved to be flawed because of the presumption of random sampling by BEAST’s FBD model. For an update on the story of this analysis, see my more recent post, Maddison’s Law of Phylogenetic Analysis.