Today, someone posted the following query to R-sig-phylo:
Does anyone know of a good implementation of an LTT plot that can draw a Confidence Interval or HPD interval from a set of trees? I've seen things like `ltt` in phytools that can draw one line for each tree in the sample. However, this can look a bit messy, and I'd ideally love to just plot the 95% CI or HPD of the ages/lineages in the trees. Has anyone seen anything like this?
Loading required package: phytools
Loading required package: paleotree
> trees<-pbtree(n=50,t=log(50)-log(2),method="direct", nsim=200,quiet=T)
This is pretty cool.
It occurred to me that there are two different CIs that we might be interested in: the CI(lineages), given a time in the past; or the CI(time) given a number of lineages. The former, CI(n|t), can be read as the (say) 95% CI of the number of lineages alive at time t; whereas the latter, CI(t|n), is the 95% CI on the timing of the nth speciation event.
Even before Dave responded, especially because phytools was mentioned in the query, I'd already started working on this. Here's a link to a function that does this, and also computes the mean (rather than median), and returns the result invisibly. So, for instance:
> # same as paleotree
time low(lineages) lineages high(lineages)
[1,] 0.00000000 1 1.0 1
[2,] 0.03218876 2 2.0 3
[3,] 0.06437752 2 2.0 3
[4,] 0.09656627 2 2.0 3
[5,] 0.12875503 2 2.0 4
[6,] 0.16094379 2 2.0 4
[7,] 0.19313255 ...
Now on time with the mode changed to "mean":
It also works for trees with varying total length or number of tips (although in the latter case, only for method="lineages". So, for instance:
Finally, we can set an arbitrary α level. For instance:
Basically, you get the idea. Please note that the version of pbtree that I'm using in the above simulations is in a non-CRAN phytools update (phytools 0.2-46).