TreeSlicer: Example 2

Setting change-point times at equally spaced intervals between the tMRCA and the present

In this example we set change-point times to be equally spaced between the tMRCA and the most recent sample (the present, as in the figure below. This is usually used for the effective reproductive number (R_e) or the birthrate. Thus, shifts in the R_e are placed at equidistant intervals between the tMRCA and the present (t₁-t₄), i.e. the time between the tMRCA (1924) and t₁ (1938) is the same as the time between each of the other shifts (e.g. t₂ and t₃, 1952-1966). Thus, only the first interval (between t₀ and t₁ is longer than the others, because it incorporates the origin branch).

Keep in mind that the time of the origin represents the start of the outbreak, which is in general older than the tMRCA of the sampled tree. (In methods based on the Kingman coalescent the origin is not estimated and the process stops at the tMRCA).

Why do this?

The default behaviour of BDSKY is to place change-point times at equally spaced intervals between the origin and the present, as below.

For most applications this is sufficient. However, when the origin branch is longer than one of the intervals, a change-point time will be placed on the origin branch (as below). This is problematic, because there is no information to identify rate-changes on the origin branch. Therefore, in some cases, equidistant change-point times between the origin and the present can cause problems with convergence and result in biased estimates.

Without TreeSlicer

It is not possible to set equidistant change-point times between the tMRCA and the present using just BDSKY, without also fixing the tree.

The closest alternative is to condition the likelihood on the tMRCA of the tree, instead of the origin of the outbreak (by setting the option conditionOnRoot).

<distribution id="BirthDeathSkySerial" 
              spec="beast.evolution.speciation.BirthDeathSkylineModel" 
              tree="@Tree"
              reproductiveNumber="@reproductiveNumber" 
              becomeUninfectiousRate="@becomeUninfectiousRate" 
              samplingProportion="@samplingProportion" 
              conditionOnRoot="true"/>

Note that in this case the origin time is not estimated.

With TreeSlicer

With TreeSlicer, we set to a TreeSlicer object, and explicitly tell it to place the change-point times between the tMRCA and the present.

<distribution id="BirthDeathSkySerial"
              spec="beast.evolution.speciation.BirthDeathSkylineModel"
              tree="@Tree"
              reproductiveNumber="@reproductiveNumber"
              becomeUninfectiousRate="@becomeUninfectiousRate"
              samplingProportion="@samplingProportion"
              origin="@origin">
    <birthRateChangeTimes spec="TreeSlicer" id="ReTreeSlice" tree="@Tree" dimension="5" stop="tmrca" includeLast="false"/>
    <reverseTimeArrays spec="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). 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.

If we had set includeLast to true in the TreeSlicer object, there would be a change-point time on the tMRCA, as below. Convergence is generally better when includeLast is set to false.

TODO: The stop and includeLast inputs need to be renamed in the next release because these names are confusing!!