TreeSlicer: Example 5

Estimating the same change-point dates across multiple trees

In this example we are interested in estimating the time of a major transition, using sequences from more than one tree as input. That is, we believe that an event independent of the trees caused a shift in one of the rates at some point in the past.

Why do this?

This may be the case when analysisng an HIV epidemic with multiple co-circulating subtypes or clusters. Implementing a single harm-reduction policy would affect all clusters. This model allows us to use data from multiple co-circulating outbreaks to estimate the time of such a major transition.

Without TreeSlicer

This may be possible without using TreeSlicer by making use of multiple RPNCalculators to transform the estimated date to a height on each of the trees. Since the trees do not generally all have the same most recent sample, a separate RPNCalculator is needed for each tree, which makes it very easy to introduce errors, especially when the data are modified in any way (since the most recent dates of each tree is now part of the model, and not restricted to the data). Because of BDSKY's idiosyncracy of requiring the change-point times to always include 0, each RPNCalculator would have to multiply by the vector (0,1) as well, in addition to performing a shift.

With TreeSlicer

Example XML: dengue4_bdsky_csc_treeslicer_joint.xml

With TreeSlicer it is possible to directly estimate the change-point time as a date, and place a prior on the date (without needing to select parts of the vector). First, add the change-point date to the state (initialised to 1970):

<parameter id="reproductiveNumberChangeDate" name="stateNode" value="1978"/>

Then change birthRateChangeTimes to a TreeDateSlicer object for both trees:

<!-- Treeprior on tree 1 -->
<distribution id="BirthDeathSkySerial.C1" 
              spec="beast.evolution.speciation.BirthDeathSkylineModel" 
              tree="@Tree.C1"
              reproductiveNumber="@reproductiveNumber.C1" 
              becomeUninfectiousRate="@becomeUninfectiousRate" 
              samplingProportion="@samplingProportion.C1" 
              origin="@origin.C1">
      <birthRateChangeTimes spec="TreeDateSlicer" id="ReTreeSlice.C1" tree="@Tree.C1" dates="@reproductiveNumberChangeDate"/>
      <reverseTimeArrays spec="beast.core.parameter.BooleanParameter" value="true false false false false"/>
</distribution>

<!-- Treeprior on tree 2 -->
<distribution id="BirthDeathSkySerial.C2" 
              spec="beast.evolution.speciation.BirthDeathSkylineModel" 
              tree="@Tree.C2"
              reproductiveNumber="@reproductiveNumber.C2" 
              becomeUninfectiousRate="@becomeUninfectiousRate" 
              samplingProportion="@samplingProportion.C2" 
              origin="@origin.C2">
      <birthRateChangeTimes spec="TreeDateSlicer" id="ReTreeSlice.C2" tree="@Tree.C2" dates="@reproductiveNumberChangeDate"/>
      <reverseTimeArrays spec="beast.core.parameter.BooleanParameter" value="true false false false false"/>
</distribution>

Now add a prior on the change-point date:

<prior id="reproductiveNumberChangeDatePrior" name="distribution" x="@reproductiveNumberChangeDate">
    <Uniform name="distr" upper="1980" lower="1975"/>
</prior>

and that's it!