Plot density of rate shifts in the entire sample

[Jamie T]

Hi,

Is it possible to extract the timings of rate shifts in every sample of the posterior for a density plot?

And the same for distinguishing whether they're up or down shifts?

[pasca...@gmail.com]

Yes, here is an example with the built-in whales dataset.

data(whales, events.whales)

ed <- getEventData(whales, events.whales, burnin=0.25, nsamples=500)

The list element ed$eventData contains a table of events for each posterior sample. This lists the timing and node for each shift, and the parameters of the exponential change model for that rate regime.

head(ed$eventData, 3)

[[1]]
  node    time     lam1       lam2        mu1 mu2 index
1   88  0.0000 0.251417 -0.0293481 0.03413370   0     1
2   42 12.8214 0.993371 -0.0464565 0.57440900   0     2
3  141 26.0480 0.301387 -0.0434084 0.00620047   0     3

[[2]]
  node    time     lam1       lam2        mu1 mu2 index
1   88  0.0000 0.180316 -0.0314398 0.00603266   0     1
2  132 16.3472 0.363055 -0.0363018 0.02545200   0     2

[[3]]
  node    time     lam1       lam2       mu1 mu2 index
1   88  0.0000 0.131678 -0.0235760 0.0103028   0     1
2  140 23.4463 0.405294 -0.0309377 0.0283852   0     2

[Jamie T]

Hi, thanks.  So based on this Dan's reply in the other thread, to determine whether a shift is an up or down shift do I need to calculate lam(t) = lam1 * exp(lam2 * t)?  I did this for a rate shift in a BSC and it comes out with a value of 2.39, but the legend in the phylorate plot spans 0.24-2.2.  Then do we need to find the direction of difference from the root regime?  How would we deal with if a rate shift is nested within a rate shift?

[pasca...@gmail.com]

The formula for lam(t) calculates the rate at any point in time within a particular rate regime, where t is the time since the start of the regime. It will not tell you whether or not that rate regime was an upward or a downward shift from the previous root-ward rate regime. BAMMtools does not have a ready-made function to do that. You would need to write some code to calculate rates just before and just after a rate shift to determine in which direction the shift went.