Significant shifts
In this vignette we will assess the support for that a branch has at least one rate shift, versus the null hypothesis of no shifts. We will test this using Bayes factors, i.e. the relative support for one or more versus no shifts.
Tree file
First, we load the necessary modules and read in the tree file.
using Pesto
phy = readtree(Pesto.path("primates.tre"))
sampling_fraction = 0.62
primates = SSEdata(phy, sampling_fraction)Model setup
Let's use the standard Pesto analysis with empirical Bayes for the speciation and extinction rates, and maximum-likelihood shift rate.
model, rates = pesto(primates; n = 6)Bayes factors
A Bayes factor is a summary of the evidence in favor of one hypothesis, as opposed to a different hypothesis (Kass & Raftery 1995). In Pesto, we use Bayes factors to assess the evidence for the hypothesis that there was at least one diversification rate shift, versus the hypothesis that there were no rate shifts, per branch. The standard equation for a Bayes factor is as follows (see also Shi & Rabosky 2015)
\[ \text{Bayes factor} = \frac{\frac{P_M(\geq \text{1 shifts})}{\pi_M(\geq \text{1 shifts})}}{\frac{P_M(0 \text{ shifts})}{\pi_M(0 \text{ shifts})}}\]
where $P_M(\geq \text{1 shifts})$ is the posterior probability of at least one shift, $P_M(\geq \text{0 shifts})$ is the posterior probability of no shifts, $\pi_M(\geq \text{1 shifts})$ is the prior probability of at least one shift, and $\pi_M(\geq \text{0 shifts})$ is the prior probability of no shifts.
Plotting Bayes factors
We can for example plot the Bayes factor directly on the tree. Since the Bayes factor can vary considerably, we instead plot the log-transformed Bayes factors, which are more concentrated around 0. A log Bayes factor with a value of 0 means that the prior and posterior support for the shift hypotheses are equal. A value much larger than 0 means that there is support for one or more shifts. A value of less than 0 means that there is more support for 0 rate shifts.
using Makie, CairoMakie
rates[!,"shift_bf_log"] = rates[!,:shift_bf]
min, max = extrema(rates[1:end-1,"shift_bf_log"])
cmap = Makie.cgrad([:gray, :black, :purple])
treeplot(primates, rates, "shift_bf_log"; cmap = cmap)
Plotting supported branches
Alternatively, we can assess which branches had a strong support for there being at least one shift. We first set the significance level at 10, meaning strong support (Kass & Raftery 1995). Next, we can compute which branches has strong support for at least one diversification rate shift, vs the null hypothesis of zero shifts.
| Bayes factor | log10 Bayes factor | Level of support |
|---|---|---|
| 0 to 3.2 | 0 to 0.5 | Not worth mentioning |
| 3.2 to 10 | 0.5 to 1 | Substantial |
| 10 to 100 | 1 to 2 | Strong |
| >100 | >2 | Decisive |
If we inspect the data frame with the branch-specific outputs, we can see specifically which branches has strong support.
using DataFrames
cutoff = 3.2
filter(:shift_bf => x -> x > cutoff, rates)| Row | delta_lambda | delta_mu | delta_netdiv | delta_relext | mean_lambda | mean_mu | mean_netdiv | mean_relext | node | edge | nshift | shift_bf | shift_bf_log |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Int64 | Int64 | Float64 | Float64 | Float64 | |
| 1 | 0.044266 | -0.03045 | 0.074716 | -0.321868 | 0.172548 | 0.0949949 | 0.0775528 | 0.580336 | 287 | 105 | 0.895757 | 7.30595 | 7.30595 |
| 2 | 0.102268 | -0.00899225 | 0.11126 | -0.236665 | 0.212724 | 0.099386 | 0.113338 | 0.519477 | 376 | 283 | 0.687401 | 5.78283 | 5.78283 |
| 3 | 0.174356 | -0.0228849 | 0.197241 | -0.400958 | 0.22018 | 0.100328 | 0.119852 | 0.548184 | 386 | 304 | 1.05207 | 39.65 | 39.65 |
The result is that there are two branches that had strong support for there being at least one shift on the branch. The branch that led to the Old World Monkeys has the highest support (Bayes factor of 6.9), while the gibbons clade (represented by Hylobates) have a slightly lower support (Bayes factor of 3.5).
In order to plot these, we can create a dummy variable and use it in the tree plotting function.
rates[!,:strong_support] = Float64.(rates[!,:shift_bf] .> cutoff)
cmap = Makie.cgrad([:black, :red])
treeplot(primates, rates, "strong_support"; cmap = cmap)
Both the phylogeny plot, as well as the filtered data frame, show two branches that had substantial support (with a Bayes factor of >3.2) for one or more shifts. The number of estimated diversification rate shifts ($\hat{N}$) on this branch is also significantly larger than zero, at 1.05 and 0.65.
References
- Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the american statistical association, 90(430), 773-795.
- Shi, J. J., & Rabosky, D. L. (2015). Speciation dynamics during the global radiation of extant bats. Evolution, 69(6), 1528-1545.