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)
```