Today a phytools user & R-sig-phylo list participant asked:
“I am performing analyses using threshBayes from the phytools package. I was wondering if someone has ever extracted the highest posterior density intervals (HPD) values using this function. How could I make this?”
As a quick reminder, the phytools function threshBayes implements a Bayesian MCMC version of Felsenstein’s (2005; 2012) threshold model for estimating the evolutionary correlation between binary character (or between a discrete & binary trait) under a phylogenetic comparative interpretation of Wright’s threshold model from quantitative genetics.
This question is a bit of a softball, since we provided an exact workflow to do precisely this in Chapter 7 of our recent Princeton University Press book; however, not everyone has the book (hard to believe, I know), so I thought it was worth reiterating here!
To illustrate this workflow I’m going to use an excellent dataset of parity mode (i.e., oviparity vs. viviparity) & environmental temperature in South American liolaemid lizards from Esquerré et al. (2019).
We can start by loading phytools.
library(phytools)
packageVersion("phytools")
## [1] '1.9.25'
Now let’s add our data & tree as follows.
data("liolaemid.data")
head(liolaemid.data)
## parity_mode max_altitude temperature
## Ctenoblepharys_adspersa O 750 23.05
## Liolaemus_abaucan O 2600 20.20
## Liolaemus_albiceps V 4020 12.38
## Liolaemus_andinus V 4900 11.40
## Liolaemus_annectens V 4688 5.10
## Liolaemus_anomalus O 1400 23.78
data("liolaemid.tree")
liolaemid.tree
##
## Phylogenetic tree with 257 tips and 256 internal nodes.
##
## Tip labels:
## Liolaemus_abaucan, Liolaemus_koslowskyi, Liolaemus_albiceps, Liolaemus_irregularis, Liolaemus_ornatus, Liolaemus_calchaqui, ...
##
## Rooted; includes branch lengths.
plotTree(liolaemid.tree,ftype="i",fsize=0.4,type="arc",
lwd=1,arc_height=0.5)