User control of arc width & color in arc.cladelabels

The other day I received the following email:

“I am playing around with the new arc.cladelabels function. Is it possible to add in options for color and line thickness. I see it calls draw.arc internally so I am assuming it would just be a matter of adding the options "col" and "lwd" in the function? It would be useful to have alternating or different colors and/or line widths for when you are labeling all clades in the tree and there is not much space between the arcs.”

Indeed, these are relatively simple options to add, which I have just done on GitHub. I also added an additional option which is to specify the node corresponding to the common ancestor of the clade to be labeled interactively.

The latter option is hard to demonstrate (absent a video); however the former options work as follows:

library(phytools)
tree

## 
## Phylogenetic tree with 200 tips and 199 internal nodes.
## 
## Tip labels:
##  t7, t23, t24, t126, t127, t35, ...
## 
## Rooted; includes branch lengths.

plotTree(tree,type="fan",lwd=1,ftype="off")
nodes

## [1] 203 215 264 288 323 356

text<-paste("Clade",LETTERS[1:length(nodes)])
text

## [1] "Clade A" "Clade B" "Clade C" "Clade D" "Clade E" "Clade F"

library(RColorBrewer)
cols<-brewer.pal(length(nodes),"Accent")
nulo<-mapply(arc.cladelabels,text=text,node=nodes,col=cols,
    MoreArgs=list(mark.node=FALSE,lwd=6))

This function version should be obtained by installing the latest development version of phytools from GitHub.

Comparing fitted discrete character models using fitMk

In the following tutorial, I'll outline how to fit & compare alternative discrete character evolution models using fitMk, a function in the phytools package that in turn utilizes code internally that is adapted from ape::ace. The same principles apply to model-fitting using geiger::fitDiscrete (which is slower, but I would recommend as more robust); however the syntax differs slightly.

First, let's see our tree & data. The tree is, of course, an object of class "phylo", and the trait is a vector of class "factor", although it could also be a character or integer vector.

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 
## c c c b c b b a c a a c b b b c c a b b b b b b b b 
## Levels: a b c

dotTree(tree,x)

Here, we'll consider six models in increasing order of parameterization: an equal-rates ("ER") model, an ordered model (but with a single transition rate), an ordered symmetric model, in which a<->b and b<->c have different rates, an ordered unsymmetric model, a symmetric, un-ordered model ("SYM"), and an all-rates-different ("ARD") model.

## ER model
fitER<-fitMk(tree,x,model="ER")
fitER

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.641937  0.320969  0.320969
## b  0.320969 -0.641937  0.320969
## c  0.320969  0.320969 -0.641937
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -24.08594 
## 
## Optimization method used was "nlminb"

plot(fitER)

## ordered single rate
model<-matrix(c(0,1,0,1,0,1,0,1,0),3,3,byrow=TRUE,
    dimnames=list(c("a","b","c"),c("a","b","c")))
fitOrdered1<-fitMk(tree,x,model=model)
fitOrdered1

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.685336  0.685336  0.000000
## b  0.685336 -1.370672  0.685336
## c  0.000000  0.685336 -0.685336
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -24.240481 
## 
## Optimization method used was "nlminb"

plot(fitOrdered1,show.zeros=FALSE)

## ordered symmetric model
model<-matrix(c(0,1,0,2,0,1,0,2,0),3,3,byrow=TRUE,
    dimnames=list(c("a","b","c"),c("a","b","c")))
fitOrdered2<-fitMk(tree,x,model=model)
fitOrdered2

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.932799  0.932799  0.000000
## b  0.688188 -1.620987  0.932799
## c  0.000000  0.688188 -0.688188
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -23.890065 
## 
## Optimization method used was "nlminb"

plot(fitOrdered2,show.zeros=FALSE)

## ordered all-rates-different model
model<-matrix(c(0,1,0,2,0,3,0,4,0),3,3,byrow=TRUE,
    dimnames=list(c("a","b","c"),c("a","b","c")))
fitOrdered4<-fitMk(tree,x,model=model)
fitOrdered4

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -1.724857  1.724857  0.000000
## b  0.543052 -1.145964  0.602912
## c  0.000000  0.811064 -0.811064
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -22.630044 
## 
## Optimization method used was "nlminb"

plot(fitOrdered4,show.zeros=FALSE)

## SYM model
fitSYM<-fitMk(tree,x,model="SYM")
fitSYM

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.406794  0.201962  0.204833
## b  0.201962 -0.669625  0.467663
## c  0.204833  0.467663 -0.672496
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -23.895651 
## 
## Optimization method used was "nlminb"

plot(fitSYM,show.zeros=FALSE)

## finally, ARD model
fitARD<-fitMk(tree,x,model="ARD")
fitARD

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           a         b         c
## a -0.710481  0.000000  0.710481
## b  0.282143 -0.829098  0.546955
## c  0.000000  0.989207 -0.989207
## 
## Fitted (or set) value of pi:
##         a         b         c 
## 0.3333333 0.3333333 0.3333333 
## 
## Log-likelihood: -22.681081 
## 
## Optimization method used was "nlminb"

plot(fitARD,show.zeros=FALSE)

Now, the question is - how can we compare among these models? Well, we can compute AICs & Akaike weights. For instance:

AIC<-setNames(sapply(list(fitER,fitOrdered1,fitOrdered2,fitOrdered4,fitSYM,
    fitARD),AIC),c("ER","Ordered-1","Ordered-2","Ordered-4","SYM","ARD"))
AIC

##        ER Ordered-1 Ordered-2 Ordered-4       SYM       ARD 
##  50.17188  50.48096  51.78013  53.26009  53.79130  57.36216

aic.w(AIC)

##         ER  Ordered-1  Ordered-2  Ordered-4        SYM        ARD 
## 0.36914685 0.31628807 0.16518555 0.07881405 0.06042989 0.01013560

This tells us that the weight of evidence is more or less evenly split between the "ER" model & our single-rate ordered model.

Alternatively, we might compare pairs of models using a likelihood-ratio test. Of course, this can only be used to compare more complex models that have simpler models as a special case. This means, for instance, that we cannot compare "ER" and "Ordered-1" because they are both equally simple & neither is a special-case of the other. We can, by contrast, compare (for example) "Ordered-1" and "Ordered-2" in this way. Let's try:

## first the number of parameters in each model:
k0<-length(fitOrdered1$rates)
k1<-length(fitOrdered2$rates)
LR<--2*(logLik(fitOrdered1)-logLik(fitOrdered2))
LR

## [1] 0.700832

P_chisq<-pchisq(LR,df=k1-k0,lower.tail=FALSE)
P_chisq

## [1] 0.4025043

That's the general idea.

Note that in this simulated case "Ordered-1" was the generating model, so we got pretty close!

Plotting a legend outside the plot area using add.color.bar (& in general)

Recently I received the following question:

“I am using the add.color.bar function to add a colorbar to an existing plot.
Do you know if it is possible to move the added colorbar to the outside of the plot domain? Do you think I would need to modify code in the actual function?”

Evidently, this can be done using par(xpd=TRUE). For instance:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

##           A           B           C           D           E           F 
##  3.09148522  1.54096027  2.91508379  2.00469068  1.15813551  0.29543463 
##           G           H           I           J           K           L 
##  0.27951056  0.60805259 -2.52405841 -3.50891295 -3.34674083 -0.97428407 
##           M           N           O           P           Q           R 
## -1.25336677 -1.40812351 -0.51371266 -2.92582374 -2.24608676 -0.53203275 
##           S           T           U           V           W           X 
## -0.59487179  0.30811419 -0.08364322  0.16895725  0.70189973 -3.64724619 
##           Y           Z 
## -3.65517946 -4.30271550

obj<-contMap(tree,x,plot=FALSE)
plot(obj,legend=FALSE,mar=c(5.1,2.1,2.1,2.1))
par(xpd=TRUE)
add.color.bar(leg=3,cols=obj$cols,title="trait value",lims=obj$lims,
    digits=3,prompt=FALSE,x=0,y=-3)

We can easily see that this is outside of our plotting domain by adding axis labels to the plot as follows:

plot(obj,legend=FALSE,mar=c(6.1,2.1,2.1,2.1))
par(xpd=TRUE)
add.color.bar(leg=3,cols=obj$cols,title="trait value",lims=obj$lims,
    digits=3,prompt=FALSE,x=0,y=-4)
axis(1)

You get the general idea.