## Thursday, August 4, 2022

### Sigmoid phylograms for "phylo", "simmap", and "contMap" object classes

Yesterday I tweeted that I'd figured out plotting a "contMap" tree (and, likewise, any object of class "simmap") in “sigmoidal” style, and promised to post more about it today.

Well, here are the functions. (I'll add this to the phytools development version shortly.)

compute.brlen.simmap<-function(phy,method="Grafen",power=1,...){
if(inherits(tree,"simmap")){
tt<-as.phylo(phy)
tt<-compute.brlen(tt) #,method=method,power=power,...)
ss<-tt\$edge.length/phy\$edge.length
phy\$maps<-mapply("*",ss,phy\$maps,SIMPLIFY=FALSE)
phy\$mapped.edge<-phy\$mapped.edge*
matrix(rep(ss,ncol(phy\$mapped.edge)),
nrow(phy\$mapped.edge),ncol(phy\$mapped.edge))
phy\$edge.length<-tt\$edge.length
} else if(inherits(tree,"phylo"))
phy<-compute.brlen(phy,method=method,power=power,...)
phy
}

sigmoidPhylogram<-function(tree,...){
## b=5,m1=0.01,m2=0.5,v=1
b<-if(hasArg(b)) list(...)\$b else 5
m1<-if(hasArg(m1)) list(...)\$m1 else 0.01
m2<-if(hasArg(m2)) list(...)\$m2 else 0.5
v<-if(hasArg(v)) list(...)\$v else 1
if(inherits(tree,"simmap")){
if(hasArg(colors)) colors<-list(...)\$colors
else {
ss<-sort(unique(c(getStates(tree,"nodes"),
getStates(tree,"tips"))))
mm<-length(ss)
colors<-setNames(
colorRampPalette(palette()[1:min(8,mm)])(mm),
ss)
}
} else if(inherits(tree,"phylo")) {
if(hasArg(color)) colors<-setNames(list(...)\$color,"1")
else colors<-setNames(par()\$fg,"1")
tree<-paintSubTree(tree,Ntip(tree)+1,"1")
}
if(hasArg(res)) res<-list(...)\$res
else res<-199
if(hasArg(use.edge.length))
use.edge.length<-list(...)\$use.edge.length
else use.edge.length<-TRUE
if(!use.edge.length){
if(hasArg(power)) power<-list(...)\$power
else power<-1
tree<-compute.brlen.simmap(tree,power=power)
}
if(hasArg(lwd)) lwd<-list(...)\$lwd
else lwd<-2
h<-max(nodeHeights(tree))
args<-list(...)
args\$power<-NULL
args\$res<-NULL
args\$colors<-NULL
args\$b<-NULL
args\$m1<-NULL
args\$m2<-NULL
args\$v<-NULL
args\$tree<-tree
args\$color<-"transparent"
dev.hold()
do.call(plotTree,args)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## Yt<-A+(K-A)/((C+exp(-B*(t-M)))^(1/v))
sigmoid<-function(x,.A=A,.K=K,.C=C,.B=B,.M=M,.v=v)
return(.A+(.K-.A)/((.C+exp(-.B*(x-.M)))^(1/.v)))
for(i in 1:nrow(tree\$edge)){
A<-pp\$yy[tree\$edge[i,1]]
K<-pp\$yy[tree\$edge[i,2]]
if(i==1) dy<-abs(A-K)
B<-b*Ntip(tree)/h
t<-seq(pp\$xx[tree\$edge[i,1]],pp\$xx[tree\$edge[i,2]],
length.out=res)
t<-sort(c(t,t[1]+cumsum(tree\$maps[[i]])))
dd<-diff(range(t))
M<-t[1] + if(m1*h>(m2*dd)) m2*dd else m1*h
C<-1
Yt<-c(A,sigmoid(t),K)
t<-c(t[1],t,t[length(t)])
COLS<-vector()
bb<-setNames(t[1]+cumsum(tree\$maps[[i]]),names(tree\$maps[[i]]))
for(j in 1:length(t))
COLS[j]<-colors[names(bb[which(t[j]<=bb)])[1]]
nn<-length(t)
segments(t[1:(nn-1)],Yt[1:(nn-1)],x1=t[2:nn],y1=Yt[2:nn],
col=COLS,lwd=lwd)
}
nulo<-dev.flush()
}

Now, let's test it out.

First, I'll graph a simple "phylo" object, in this case a tree of primates from Kirk & Kay (2004).

library(phytools)
primate.tree<-read.tree(file="http://www.phytools.org/Rbook/3/primateEyes.phy")
sigmoidPhylogram(primate.tree,ftype="i",fsize=0.5,lwd=1)

Next, let's see a stochastic mapped tree.

For this, I'll pull the corresponding primate data set for activity pattern, skull size, and eye size.

primate.data<-read.csv(file="http://www.phytools.org/Rbook/3/primateEyes.csv",
row.names=1,stringsAsFactors=TRUE)
head(primate.data)
##                                  Group Skull_length Optic_foramen_area Orbit_area
## Allenopithecus_nigroviridis Anthropoid         98.5                7.0      298.7
## Alouatta_palliata           Anthropoid        109.8                5.3      382.3
## Alouatta_seniculus          Anthropoid        108.0                8.0      359.4
## Aotus_trivirgatus           Anthropoid         60.5                3.1      297.4
## Arctocebus_aureus            Prosimian         49.5                1.2      134.8
## Arctocebus_calabarensis      Prosimian         53.8                1.6      156.6
##                             Activity_pattern Activity_pattern_code
## Allenopithecus_nigroviridis          Diurnal                     0
## Alouatta_palliata                    Diurnal                     0
## Alouatta_seniculus                   Diurnal                     0
## Aotus_trivirgatus                  Nocturnal                     2
## Arctocebus_aureus                  Nocturnal                     2
## Arctocebus_calabarensis            Nocturnal                     2
Activity_pattern<-setNames(primate.data\$Activity_pattern,
rownames(primate.data))
primate.maps<-make.simmap(primate.tree,Activity_pattern,model="ARD",
pi="fitzjohn",nsim=100)
## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
##              Cathemeral      Diurnal    Nocturnal
## Cathemeral -0.050529521  0.050529521  0.000000000
## Diurnal     0.002306398 -0.004367451  0.002061052
## Nocturnal   0.002586387  0.003565036 -0.006151422
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##  Cathemeral     Diurnal   Nocturnal
## 0.003675779 0.021703763 0.974620458
## Done.
cols<-setNames(c("#87CEEB","#FAC358","black"),levels(Activity_pattern))
sigmoidPhylogram(primate.maps[[1]],colors=cols,ftype="i",fsize=0.5)

It's actually no trouble for us to overlay node probabilities from our full set of stochastic maps onto this tree.

Here, I'll do it but plot only the nodes with ambiguity (posterior probability < 0.95) for the node state.

sigmoidPhylogram(primate.maps[[1]],colors=cols,ftype="i",fsize=0.5)
node.probs<-summary(primate.maps)\$ace[1:primate.tree\$Nnode,]
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
h<-max(nodeHeights(primate.tree))
nn<-Ntip(primate.tree)
for(i in 1:primate.tree\$Nnode)
if(max(node.probs[i,])<0.95)
plotrix::floating.pie(pp\$xx[i+nn],pp\$yy[i+nn],
x=node.probs[i,],radius=0.02*h,col=cols,
border=FALSE)
legend("topleft",names(cols),pch=15,col=cols,pt.cex=1.5,
cex=0.8,bty="n")

So far, so good.

Lastly, as I tweeted yesterday, we also apply this to a "contMap" style graph.

To do this, we'll need to simply extract the tree from the "contMap" object class (as this is a "simmap" object) and pass it to sigmoidPhylogram.

To illustrate this, I'll use a phylogeny & data for cordylid lizards from Broeckhoven et al. (2016).

In this example, I'm going to plot my phylogeny twice – to create the effect of an outline – and I'll add a color legend at the end.

Readers will see that I also adjusted one argument: m1. This (along with b, m2, and v) affect how the sigmoids are plotted, but have already been optimized so that the tree (generally) plots pretty well with the defaults. Try different values and see what happens!

cordylid.tree<-read.tree(file="http://www.phytools.org/Rbook/5/cordylid-tree.phy")
cordylid.data<-read.csv(file="http://www.phytools.org/Rbook/5/cordylid-data.csv",
row.names=1)
cordylid.pc1<-setNames(cordylid.data\$pPC1,rownames(cordylid.data))
cordylid.cMap<-contMap(cordylid.tree,cordylid.pc1,plot=FALSE)
cordylid.cMap<-setMap(cordylid.cMap,terrain.colors(n=10))
sigmoidPhylogram(cordylid.tree,ftype="i",fsize=0.9,b=5,m1=0.05,
lwd=7)
sigmoidPhylogram(cordylid.cMap\$tree,colors=cordylid.cMap\$cols,
ftype="off",b=5,m1=0.05,lwd=5,add=TRUE,
xlim=get("last_plot.phylo",envir=.PlotPhyloEnv)\$x.lim,
ylim=get("last_plot.phylo",envir=.PlotPhyloEnv)\$y.lim)
add.color.bar(0.5*max(nodeHeights(cordylid.tree)),
cols=cordylid.cMap\$cols,title="phylogenetic PC 1",
lims=cordylid.cMap\$lims,digits=2,prompt=FALSE,lwd=5,
x=0,y=1,subtitle="",fsize=0.8)