Monday, October 30, 2023

Computing an HPD interval on the evolutionary correlation from threshBayes in phytools

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)