Allow calibration of tips

[Anaphory]

While my criticism of the calibration in #89 was hasty, I think one of its points remains.

We currently do not allow calibration of a tip, but I think that's a very frequent use case for calibrations. This issue is actually two-fold (maybe three, I need to test a beast2 thing.)

1. If a given clade matches only one language, a test of if len(langs) < 2: will still result in

   [INFO] Calibration on clade macr1272 ignored as no matching languages in analysis.

2. get_languages_by_glottolog_clade(clade) does not actually step to the tips, which is why the previous example cites macr1272, which I used when fishing for port1283 only – Specifying port1283 as the clade to calibrate should work, but leads to get_languages_by_glottolog_clade(clade) returning [].
3. (Possibly) I think specifying tip dates in Beast may not be as easy as “prior on this single-taxon TaxonSet's MRCA”, but it may be.

[Anaphory]

I'm really not good at thinking straight today, am I? First #89, then this. Of course, there is an even bigger issue here: There actually needs to be an operator changing the tip-height for this to work, and to start with it would be useful to have the restricted tip starting at an appropriate age, so the <plate> around the taxon definitions would have to change (or the tip date set in some other fashion).

[Anaphory]

Some changes to do this now in my fork, in tip_calibrations: https://github.com/Anaphory/BEASTling/commit/69d473000431edf175a47edd0d42cb00d3cf3dfd

[Anaphory]

Actually, I had not noticed this before:

WARNING: Yule Model cannot handle dated tips. Use for example a coalescent prior instead.

Maybe it's time to re-visit how we build the tree…

(Also, this throws a spanner in my idea of comparing different models with reversible jumps on a “trusted” dated tree to find out how languages really behave.)

[lmaurits]

Am I right in extracting the take home message that there's nothing to be done about this issue for as long as BEASTling only supports the Yule prior? I do think eventually we will have to migrate away from that restriction, but not for the next release.

[Anaphory]

Sounds good to me.

[lmaurits]

I promised @LinguList I would fast-track this feature, as it's the only thing missing from BEASTling that he needs to run a Sino-Tibetan analysis with data for an extinct Burmese language. I am keen to support papers using BEASTling by authors who aren't core devs, so it would be good to see this in the next release.

I don't think this is hard at all. The tree prior is generated in a very well-defined location, all we need is to set a flag when using tip calibrations that trigers the use of the coalescent prior. Perhaps we could expose this directly via a tree_prior option in [languages], but I'm not so sure about that. At the very least if somebody specifies a coalescent prior and they are not using tip calibrations, BEASTling should ask them very loudly if they know what they are doing.

I would say leave that ability out entirely, but I actually do have ideas for other tree priors I want to support in the future (but not as soon as 1.4), so that option is going to appear eventually anyway.

[Anaphory]

Yes, please make the prior choice explicit. I think it would be strange to get an entirely different tree prior depending on whether all languages are contemporary or one of them is from a year ago and all others are contemporary. (This is in addition to the later option of adding other tree priors.)

And if tree priors implicitly switch depending on some other feature, the documentation should be very loud about that.

[lmaurits]

Well, I had envisaged the switch from Yule to coalescent in the case of calibrated tips being implicit (although with a --verbose notification, for sure), because it is basically required for the analysis to be mathematically valid (for now). I don't want it to be easy (or perhaps even possible) to create invalid analyses with BEASTling, but I would also prefer that achieving this doesn't require users to know lots of exotic things about tree priors (or, in another case, about ascertainment correction), so implicintness seems somewhat unavoidable...

Regarding the ability to deliberately create things with BEASTling that are "wrong": I've wanted to do this a few times myself lately just for the sake of testing things. I've been toying with adding a god_mode = True option to [admin] which enables some other settings that I would otherwise not want people to touch, but perhaps this is not a good idea and I should just maintain a private branch in my own working directory where I can add crazy options.

[lmaurits]

@Anaphory I have just merged the tip_calibration branch from your fork into the develop branch of the main repo, so future work on this should happen there.

(I may take care of the rest of this later today, if I have time)

[lmaurits]

Notes to myself:

• If you set a tip height for just one taxon in a TraitSet because you want all other taxa to have the standard height of zero, BEAST seems to actually assume that all taxa have the height you provided for your one extinct language. So I think we need to explicitly set all others to zero?! Look into this...
• There appears to be two operators for changing tip heights, TipDatesScaler and TipDatesRandomWalker. I have not managed to get the Scaler to work at all (the tip date never changes from the initial value, regardless of parameter settings), whereas the RandomWalker appears to work just fine (and is the only operator mentioned in the official tutorial).
• As far as I can see, the sampled height of a tip does not actually affect the logged Newick tree at all, i.e. ultrametric trees still come out the other end. I guess the branch length leading to taxa with sampled heights are just truncated for the sake of likelihood calculations. It's possible we need to let the Newick logger know about the dates TraitSet in order for it to output the correct Newicks - or not, I need to look into this. For now, the only way to be sure tip dates are being sampled and to see the estimated height is to log the tip date.

[lmaurits]

Another note: The BEAST FAQ has an entry on tips sampling and seems to suggest one operator is required per calibration. I can't understand why, even after looking at the operator code. Starting to wonder if thta approach is just because of the assumption that you are trying to make minimal hand edits to a BEAUTi-produced file, and copying and pasting an operator and adding an index is easier than defining a custom TaxonSet?

[lmaurits]

Ah, it's because it moves all taxa in the taxonset in lockstep, rather than moving a random one each time. So we really do need one per calibration. The implementation I just pushed will have to be changed, then. Next week!

[lmaurits]

Specifying a tip calibration now causes a constant-population coalescent prior to be used in place of a Yule prior. This runs, but in the toy examples I've played with to ensure it runs the population size parameter invariably blows up. Not sure if this is a genuine problem or just happens because the data in tests/data/basic.csv probably carries very weak signal. Need to look into this a little more.

There is now one TipDatesRandomWalker per scaled tip, as is required. There seems to be no problem with multiple tips calibrated to different date ranges getting correctly sampled.

Trying to get to the bottom of the issue mentioned above where taxa without an explicit trait setting get initialised to some non-zero value has lead me on a magical journey into the heart of BEAST to see how this tip date stuff is actually implemented. I have found some answers and some questions.

In the way of answers, the date TraitSet can be specified as being of type date, date-forward or date-backward. The first two seem to function identically. As far as I can tell, date-forward is intended for actual literal dates (like 2017, 1985, etc.). In this case, taxa with no explicit date are set to the maximum value of the explicit dates, i.e. if three taxa were given the dates 2017, 1985 and 1970 then BEAST would assume that all of the other taxa are from 2017. I think date-backward is for expressing times in years BP (or any unit BP), and in this case taxa with no implicit date are initialised to zero (i.e. the present). I think this is the behaviour we have in mind for BEASTling (as non-tip calibrations work in exactly this way), so probably we need to change to specifying our date TraitSet as a date-backward type.

In the way of questions is how does any of this end up affecting the likelihood calculation? I've learned that BEAST does not explicitly store branch lengths. Rather, each node of the tree has a "height" stored, which measures the sum of branch lengths to a descendant child. Leaf nodes typically have this set to zero. When computing likelihoods, branch lengths are computed on the fly by subtracting the parent's height from the child's height. The TipDatesRandomWalker operator just alters a tip's height value. This all explains why the TreeLikelihood classes don't need to know anything about the date TraitSet. Presumably when the tree is initialised, something sets all the leaf node heights to be consistent with the provided date traits. I just haven't yet figured out where that happens. Once I do I should have some clearer understanding of how the different date types are interpreted....

[lmaurits]

My earlier report that BEAST logs ultrametric trees when sampling tip dates appears to have been some kind of mistake/confusion. Perhaps I was looking at trees sampled before I figured out that the TipDatesScaler wasn't doing anything. So, there's no problem there.

[Anaphory]

Ah! I think I tried to figure that out once, as well, with the added complication that I tried it while translating a file from Beast1 to Beast2, which does not make it any easier. I cannot remember whether I gained the necessary enlightenment, in any case I don't have it any more.

[lmaurits]

(assuming you were talking about initialising node heights) I figured it out in the end, the logic is a little spread out between the Tree class and the Node class, whereas I was only looking at Tree when I wrote the above. I am thinking of starting a repo of Markdown text files documenting these kinds of deep dives into undocumented BEAST internals, as a kind of public service.

I am now convinced that using date-backwards is required for BEAST to interpret BEASTling tip calibrations in the way we are intending to interpret them and have changed the code accordingly. Taxa without explicit dates are now automatically set to zero rather than the highest calibration, as expected.

There is a question around whether we want to also support the other kind of date formatting (i.e. providing years). While it bothers me a little that having to specify dates in units BP means that config files eventually "get old" (although this is probably not a practical problem, over the span of a typical academic career), I am leaning pretty strongly toward not supporting absolute dating because it will only lead to things like AD vs BC headaches.

[lmaurits]

Recent commit 714442f78b20e95c12c9ba073d9de874481d8d93 made a change to the logic surrounding how we decide which calibrations should be interpreted as tip calibrations, because I think the previous approach was possibly too lenient. The logic behind the new code should hopefully be quite clear from a quick read, I would appreciate a quick LGTM from anybody who has time.

[lmaurits]

Re: earlier reports of exploding populationSize parameters in analyses with tip calibrations. I am now running an analysis using real data, with some MRCA calibrations in addition to the tip calibration. I am not seeing any parameter explosions, the popSize seems to be holding steady within a finite region, the tree height is converging to where I expect it, and the tip is slowly moving toward the centre of its calibrated range. I think the coalescent prior implementation is probably fine, and the problems I saw earlier were caused by a combination of very small taxon count, synthetic data with no strong phylogenetic signal and extreme lack of timing information coming from only one calibration in the form of a wide interval on the height of a single tip.

[lmaurits]

I just did another real data test of this using two tip calibrations with non-overlapping ranges and everything seemed to work fine. All of the problems I mentioned encountering along the way above are now solved as far as I can tell, and I just committed documentation for this, so I think that - after over a year in waiting! - this can be closed now.