TreeSlicer: Example 4

Estimating change-point times of rates as dates

In this example we are interested in estimating the time of a major transition. That is, we are estimating a single change-point time between the origin of the tree and the present.

Why do this?

This is useful for testing hypotheses about times of major transitions, such as:

• Did the reproductive number change after the introduction of harm-reduction measures?
• Is the infectious period shorter after the introduction of treatment?
• Did the speciation rate increase after the K-T extinction?

Without TreeSlicer

It is possible to estimate the time of a single parameter shift using BDSKY. First, add the change-point times to the state:

<parameter id="reproductiveNumberChangeTime" name="stateNode" dimension="1">0 24</parameter>

Note that this vector needs to be of length 2 for estimating 1 change-point time (BDSKY requires a 0 in the change-point times vector for a reason I don't quite understand). We initialise the change-point time to 1970 (the most recent sampling time is 1994, 1994-1970=24).

Next, we add this parameter to the birth-death skyline tree-prior as the birthRateChangeTimes:

<distribution id="BirthDeathSkySerial" 
              spec="beast.evolution.speciation.BirthDeathSkylineModel" 
              tree="@Tree"
              reproductiveNumber="@reproductiveNumber" 
              becomeUninfectiousRate="@becomeUninfectiousRate" 
              samplingProportion="@samplingProportion" 
              origin="@origin">
    <birthRateChangeTimes idref="reproductiveNumberChangeTime"/>
    <reverseTimeArrays spec="beast.core.parameter.BooleanParameter" value="true false false false false false"/>
</distribution>

We also need to indicate to BDSKY that the reproductive number (or birth-rate) change-point times should be reversed, meaning that change-point times are measured as time between the present and the change-point time (heights on the tree), instead of between the origin and the change-point time (default behaviour for BDSKY). Measuring the change-point time from the origin would lead to possibly poor convergence, because the change-point time will change every time the origin time is changed. We do this using the reverseTimeArrays array. Reproductive number (or birth-rate) change-point times are the first in the array (see BDSKY code/documentation), which is why the first element is set to true.

To set a prior on the change-point time we first have to select only the second element in the change-point times vector and then place a prior on it alone. This is because the prior probability is evaluated at every element of a vector, and we don't necessarily want to include 0 in our prior on the change-point height (meaning that the present, 1994 in this case, is included in the prior for the change-point date). To do this we can use beast.core.util.Slice or ExcludablePrior (can be found in multitypetree.distributions or beast.math.distributions). I prefer to use Slice, since it is part of the BEAST2 core and it involves less typing. First declare the slice somewhere in the XML file (not inside the MCMC object):

<function spec="beast.core.util.Slice" id="reproductiveNumberChangeTimeSlice" arg="@reproductiveNumberChangeTime" index="1" count="1"/> 

This tells BEAST2 to select 1 element from reproductiveNumberChangeTime, starting at element 1 (Slice starts number at 0). Now we can place a prior on the change-point time, in this case we add a uniform prior between 1960 and 1980. Since the most recent sample is from 1994, the uniform prior is between 14 and 34 (note that the upper and lower limits reverse when switching from dates to heights):

<prior id="reproductiveNumberChangeTimePrior" name="distribution" x="@reproductiveNumberChangeTimeSlice">
    <Uniform name="distr" upper="34" lower="14"/>
</prior>

With TreeSlicer

Why use TreeSlicer?

So why use TreeSlicer?