Skip to content

TreeSlicer: Example 6

Louis edited this page Aug 27, 2019 · 2 revisions

Setting the same change-point dates across multiple trees

In this example we set the change-point times for two trees to be the at the same dates. We set the change-point times of the first tree to be equally spaced between the tMRCA and the most recent sample. The change-point times for the second tree are then calculated to fall at the same dates. Because change-point times cannot be negative, we choose the first tree to be the tree with the oldest most-recent sample.

Why do this?

When there are multiple co-circulating clusters of an infectious disease, it may be that some epidemiological parameters are influenced by all clusters (e.g. the reproductive number, through density dependence). In this case a joint analysis using data from all clusters would be more informative than independent analyses. Using TreeSlicer to use the same change-point dates across multiple trees, without conditioning on the dates, allows us to do this efficiently.

Without TreeSlicer

I don't know if it is possible to set up a model like this using just BDSKY alone. I suspect not, but if it is, it would be a mess of RPNCalculators and Slices. It would be incredibly easy to inadvertently introduce errors into such an analysis.

With TreeSlicer

Example XML: dengue4_bdsky_equidistant5_treeslicer_joint.xml

With TreeSlicer, we change birthRateChangeTimes of the first tree (the tree with the oldest most-recent sample) to a TreeSlicer object, and explicitly tell it to place the change-point times between the tMRCA and the present.

<!-- Treeprior of 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="TreeSlicer" id="ReTreeSlice.C1" tree="@Tree.C1" dimension="3" to="tmrca" inclusive="false"/>
                        <reverseTimeArrays spec="beast.core.parameter.BooleanParameter" value="true false false false false"/>
                  </distribution>

We add a TreeSliceDateLogger to the tracelog, to transform the change-point times to dates (and log them for post-processing, logging the dates to a file is not necessary, but is generally good practice and makes it easier to reconstruct the marginal posterior distribution at different timepoints).

<log id="ReSliceDates.C1" spec="TreeSliceDateLogger" treeSlice="@ReTreeSlice.C1"/>

Next, we need to use these dates as input to the change-point times for the second tree. However, these dates include the most-recent sampling date of tree 1, which we want to remove. We use a Slice object to slice it off.

<function spec="Slice" id="ReSliceDates.C1.trim" arg="@ReSliceDates.C1" index="1" count="2"/>

Now we can specify the change-point times of the second tree using a TreeDateSlicer, and use the output of the above Slice, which contains the change-point dates of the first tree, as input dates.

<!-- Treeprior of 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="@ReSliceDates.C1.trim"/>
                        <reverseTimeArrays spec="beast.core.parameter.BooleanParameter" value="true false false false false"/>
                  </distribution>

We can also log the change-point dates of the second tree as a sanity check, to make sure that both trees have the same change-point dates at every step of the MCMC algorithm. (Note: the most recent date logged for each tree will be the most recent sampling date of each tree, and not the first change-point time, and so will be different between the two trees).

<log id="ReSliceDates.C2" spec="TreeSliceDateLogger" treeSlice="@ReTreeSlice.C2"/>