Generic simmap method for "simmap" object class. (Wait... what?)

As documented on this blog, and in my phytools 2.0 pre-print, I recently added a new generic simmap method for stochastic mapping that wraps the traditional phytools::make.simmap function. The idea of this generic is to make it easier for phytools users to (for instance) fit multiple discrete character evolution models using fitMk, and then run stochastic mapping on the best of these (without having to re-estimate Q), or even to sampling across alternative models in proportion to their model weights.

Today I have also added a generic simmap method for the "simmap" object class.

Wait – what?

The way this method is designed is that it pulls the tip states (and Q transition matrix, if available) out of an existing "simmap" object, and then runs stochastic mapping with these data.

Here’s how a simple demo using a phylogeny of tropidurid lizards with a mapped character (rock vs. non-rock dwelling), as used in Revell et al. (2022).

## load phytools
library(phytools)
packageVersion("phytools")

## [1] '1.5.15'

## load a "simmap" tree
data("tropidurid.tree")
tropidurid.tree

## 
## Phylogenetic tree with 76 tips and 75 internal nodes.
## 
## Tip labels:
## 	Leiocephalus_raviceps, Leiocephalus_carinatus, Leiocephalus_psammodromus, Leiocephalus_personatus, Leiocephalus_barahonensis, Stenocercus_ochoai, ...
## 
## The tree includes a mapped, 2-state discrete character
## with states:
## 	n_rock, rock
## 
## Rooted; includes branch lengths.

Let’s plot this tree.

cols<-setNames(c("white","black"),c("n_rock","rock"))
sigmoidPhylogram(tropidurid.tree,direction="upwards",fsize=0.5,
  ftype="i",colors=cols,lwd=3,outline=TRUE)
legend("bottomright",c("non-rock","rock"),pch=22,pt.bg=cols,
  pt.cex=1.5,cex=0.9,bty="n")

Great, now let’s run our generic method on this phylogeny, e.g.:

trop_trees<-simmap(tropidurid.tree,model="ARD")
trop_trees

## 100 phylogenetic trees with mapped discrete characters

trop_dmap<-densityMap(trop_trees,plot=FALSE)

## sorry - this might take a while; please be patient

trop_dmap<-setMap(trop_dmap,cols)
plot(trop_dmap,direction="upwards",fsize=c(0.5,0.9),
  ftype="i",colors=cols,lwd=3,outline=TRUE,offset=1,
  legend=1)

It works!