Additional user control for plotting "multiSimmap" densities

Today, a phytools user asked:

PCM/evolution twitter: im looking to change the colors on this density plot i made from my simmaps. i tried the base R way but it doesn't seem to work (col = c(col1, col2)). halp pic.twitter.com/REhj4jjeOe

— kent sorgon (rage era) (@kntsrgn) May 18, 2022

Indeed, it's true that the plot method for objects of this class is not particularly flexible. At the time that I wrote it, I was more focused on making something that looked good & worked, rather than code that passed maximum control to the user.

Kent's problem, however, is easy enough to fix, so I just pushed a phytools update to GitHub with the suggested improvement.

Here's a quick example using a phylogeny and dataset for presence/absence of spines from Ramm et al. (2019).

First, some preliminaries.

## load phytools and check version number
library(phytools)
packageVersion("phytools")

## [1] '1.0.8'

## read tree and data from file
lizard.tree<-read.nexus(file=
    "http://www.phytools.org/Rbook/7/lizard_tree.nex")
lizard.tree

## 
## Phylogenetic tree with 4162 tips and 4161 internal nodes.
## 
## Tip labels:
##   Sphenodon_punctatus, Dibamus_bourreti, Dibamus_greeri, Dibamus_montanus, Anelytropsis_papillosus, Dibamus_tiomanensis, ...
## 
## Rooted; includes branch lengths.

## read data from file
lizard.data<-read.csv(file=
    "http://www.phytools.org/Rbook/7/lizard_spines.csv",
    stringsAsFactors=TRUE,row.names=1)
head(lizard.data)

##                                   habitat tail.spines
## Ablepharus_budaki          non-saxicolous   non-spiny
## Abronia_anzuetoi           non-saxicolous   non-spiny
## Abronia_frosti             non-saxicolous   non-spiny
## Acanthodactylus_aureus     non-saxicolous   non-spiny
## Acanthodactylus_opheodurus non-saxicolous   non-spiny
## Acanthodactylus_pardalis   non-saxicolous   non-spiny

## load geiger and prune tree to include only species
## with data
library(geiger)
chk<-name.check(lizard.tree,lizard.data)
summary(chk)

## 3504 taxa are present in the tree but not the data:
##     Ablepharus_chernovi,
##     Ablepharus_kitaibelii,
##     Ablepharus_pannonicus,
##     Abronia_aurita,
##     Abronia_campbelli,
##     Abronia_chiszari,
##     ....
## 
## To see complete list of mis-matched taxa, print object.

lizard.tree<-drop.tip(lizard.tree,chk$tree_not_data)
## extract discrete character
tail.spines<-setNames(lizard.data$tail.spines,
    rownames(lizard.data))

Now my stochastic character mapping.

## stochastic mapping
simmap<-make.simmap(lizard.tree,tail.spines,nsim=100,
    model="ARD",pi="fitzjohn")

## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##              non-spiny        spiny
## non-spiny -0.001235890  0.001235890
## spiny      0.003977541 -0.003977541
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##  non-spiny      spiny 
## 0.98834922 0.01165078

## Done.

simmap

## 100 phylogenetic trees with mapped discrete characters

## generate a plot
dmap<-densityMap(simmap,plot=FALSE)

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

dmap<-setMap(dmap,c("yellow","darkred"))
par(bg="black",fg="white")
plot(dmap,lwd=c(1,6),ftype="off")

Now let's compute our "changesMap" object and plot it using our new control of color. Here, I decided to have the user separately set the colors and the transparency level, alpha. alpha = 0 corresponds to full transparency (i.e., invisible!), while alpha = 1 should be opaque.

## compute density
density<-density(simmap)
## graph density
par(las=1)
plot(density,colors=c("darkred","yellow"),alpha=0.5)
grid()
box()

Awesome. That's totally what we were going for!