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.

OK, I can't believe that worked.#Rstats #phytools continuous character (“contMap”) plot graphed in a sigmoidal phylogram style.

Code tomorrow. pic.twitter.com/gtqdywuihA

— Liam Revell (@phytools_liam) August 4, 2022

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)

Come on. That's not too bad, right?