Slides for Evolution 2017: A phytools plotting demo

The following is a demo of the phytools functions & methods I used to create the graphics for my Evolution 2017 meeting talk entitled “Some additional new methods for plotting phylogenies & comparative data.” The final slides for this talk are available here.

library(phytools)
packageVersion("phytools")

## [1] '0.6.18'

## -- Read some data from file:

tree<-read.tree("tree.tre")
x<-as.matrix(read.csv("x.csv",row.names=1))[,1]

## -- These methods simply plot the data at the tips of the tree. There also
## -- exists a number of other comparable functions.

dotTree(tree,x,length=10,ftype="i")

plotTree.barplot(tree,x)

plotTree.barplot(tree,x,args.barplot=list(col=sapply(x,
    function(x) if(x>=0) "blue" else "red"),
    xlim=c(-4,4)))

## -- We can also plot the data at the tips of the tree for multiple traits
## -- simultaneously. For instance:

X<-read.csv("X.csv",row.names=1)
dotTree(tree,X,standardize=TRUE,length=10)

phylo.heatmap(tree,X,standardize=TRUE)

## -- Aside from plotting the data at the tips, we can also visualize its
## -- reconstructed evolution. A simple method for doing this is called a 
## -- 'traitgram' in which the phylogeny is projected into phenotype space
## -- as follows:

phenogram(tree,x,spread.labels=TRUE,spread.cost=c(1,0))

## -- A variant of this function exists where the idea is to visualize the 
## -- uncertainty around the reconstructed values at internal nodes, as
## -- follows:

fancyTree(tree,type="phenogram95",x=x,spread.cost=c(1,0))

## Computing density traitgram...

## -- Another method in phytools projects reconstructed states onto internal
## -- edges & nodes of the tree using a color gradient. I will also demonstrate
## -- a helper function called errorbar.contMap which computes & adds error
## -- bars to the tree using the same color gradient as was employed in the 
## -- plot.

obj<-contMap(tree,x,plot=FALSE)
plot(obj,lwd=7,xlim=c(-0.2,3.6))
errorbar.contMap(obj)

## -- This method creates a special object class that can then be replotted
## -- without needing to repeat all the computation. We can also replot this
## -- object using different parameters. For example, in the following script
## -- I will invert & then gray-scale the color gradient of our original plot
## -- & then replot these two trees in a facing manner.

obj

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (-4.964627, 3.299418).

obj<-setMap(obj,invert=TRUE)
grey<-setMap(obj,c("white","black"))
par(mfrow=c(1,2))
plot(obj,lwd=7,ftype="off",xlim=c(-0.2,3.6),legend=2)
plot(grey,lwd=7,direction="leftwards",ftype="off",xlim=c(-0.2,3.6),
    legend=2)

## -- Now let's try the same thing but with some empirical data:

anole.tree<-read.tree("Anolis.tre")
anole.tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
## 
## Rooted; includes branch lengths.

svl<-read.csv("svl.csv",header=TRUE,row.names=1)
svl<-setNames(svl[,1],rownames(svl))
obj<-contMap(anole.tree,svl,plot=FALSE)
obj<-setMap(obj,invert=TRUE)
plot(obj,fsize=c(0.4,1),outline=FALSE,lwd=c(3,7),leg.txt="log(SVL)")

## -- This plotting method is fairly flexible, as we've seen. We can also
## -- plot our tree in a circular or "fan" style. In the following demo,
## -- I will also show how to add raster images of some of the taxa to
## -- our plot, which is kind of cool.

plot(obj,type="fan",outline=FALSE,fsize=c(0.6,1),legend=1,
    lwd=c(4,8),xlim=c(-1.8,1.8),leg.txt="log(SVL)")
spp<-c("barbatus","porcatus","cristatellus","equestris","sagrei")
library(png)
w<-0.5
for(i in 1:length(spp)){
    xy<-add.arrow(obj,spp[i],col="transparent",arrl=0.45,lwd=3,hedl=0.1)
    arrl<-if(spp[i]%in%c("sagrei","porcatus")) 0.24 
        else if(spp[i]=="barbatus") 0.2 else 0.34
    add.arrow(obj,spp[i],col="blue",arrl=arrl,lwd=3,hedl=0.05)
    img<-readPNG(source=paste(spp[i],".png",sep=""))
    asp<-dim(img)[1]/dim(img)[2]
    rasterImage(img,xy$x0-w/2,xy$y0-w/2*asp,xy$x0+w/2,xy$y0+w/2*asp)
}

## -- Another popular visualization method in phytools is something called a
## -- phylomorphospace, which is just a projection of the tree into (normally)
## -- a two-dimensional morphospace. For example:

phylomorphospace(tree,X[,1:2],xlab="trait 1",ylab="trait 2")

## -- For more than three continuous trait dimensions phytools also has a method
## -- that I really like, although I've literally never seen it used in 
## -- publication. I call this method a 'phylogenetic scatterplot matrix' 
## -- because it contains unidimensional contMap plots on the diagonal &
## -- bivariate phylomorphospaces for each pair of traits in off-diagonal 
## -- positions. This is what that looks like for five traits:

obj<-fancyTree(tree,type="scattergram",X=X[,1:5])

## Computing multidimensional phylogenetic scatterplot matrix...

## -- In addition to these continuous character methods, phytools also has a
## -- number of discrete character visualization tools for plotting discrete
## -- valued traits & trees.

## -- For instance, phytools contains a range of different methods related
## -- to a statistical approach called stochastic character mapping, in which
## -- dicrete character histories are sampled from their posterior probability 
## -- distribution. For starters, I will just use a stochastic character map
## -- tree that can be loaded from the package.

data(anoletree)
cols<-setNames(palette()[1:6],mapped.states(anoletree))
plot(anoletree,cols,type="fan",fsize=0.8,lwd=3,ftype="i")
add.simmap.legend(colors=cols,x=0.9*par()$usr[1],
    y=0.9*par()$usr[4],prompt=FALSE,fsize=0.9)

## -- Normally, we would be working with a set of stochastic character
## -- histories, not just a single such history. Here's an example in which
## -- we first do stochastic mapping, & then we plot the posterior 
## -- probabilities that each node is in each state.

ecomorph<-as.factor(getStates(anoletree,"tips"))
trees<-make.simmap(anoletree,ecomorph,model="ER",nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##             CG          GB          TC          TG          Tr          Tw
## CG -0.11570723  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145
## GB  0.02314145 -0.11570723  0.02314145  0.02314145  0.02314145  0.02314145
## TC  0.02314145  0.02314145 -0.11570723  0.02314145  0.02314145  0.02314145
## TG  0.02314145  0.02314145  0.02314145 -0.11570723  0.02314145  0.02314145
## Tr  0.02314145  0.02314145  0.02314145  0.02314145 -0.11570723  0.02314145
## Tw  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145 -0.11570723
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##        CG        GB        TC        TG        Tr        Tw 
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

## Done.

obj<-summary(trees,plot=FALSE)
plot(obj,colors=cols,type="fan",fsize=0.8,cex=c(0.5,0.3))
add.simmap.legend(colors=cols,x=0.9*par()$usr[1],
    y=0.9*par()$usr[4],prompt=FALSE,fsize=0.9)

## -- These are not visualizations for phylogenies, per se, but there are 
## -- several other related plotting methods in the phytools that we can 
## -- employ here. For instance, we can visualize our fitted discrete 
## -- character model:

fit.ARD<-fitMk(anoletree,ecomorph,model="ARD")
plot(fit.ARD,main="fitted ARD model for ecomorph evolution")

plot(fit.ARD,show.zeros=FALSE,main="fitted ARD model for ecomorph evolution")

## -- We can also visualize the posterior density of changes of different types
## -- on the tree. For example:

dd<-density(trees)
plot(dd)

## -- For binary traits we have more options still. In this example we will
## -- use a different empirical dataset consisting of feeding mode (biting vs.
## -- suction feeding) in elomorph eels. The first plot is just the character
## -- distribution at the tips; whereas the second is the posterior density
## -- from stochastic mapping at the edges & nodes of the tree.

eel.tree<-read.tree("elopomorph.tre")
eel.data<-read.csv("elopomorph.csv",row.names=1)
fmode<-as.factor(setNames(eel.data[,1],rownames(eel.data)))
dotTree(eel.tree,fmode,colors=setNames(c("blue","red"),
    c("suction","bite")),ftype="i",fsize=0.7)

eel.trees<-make.simmap(eel.tree,fmode,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##                bite     suction
## bite    -0.01582783  0.01582783
## suction  0.01582783 -0.01582783
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##    bite suction 
##     0.5     0.5

## Done.

obj<-densityMap(eel.trees,states=c("suction","bite"),plot=FALSE)

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

plot(obj,lwd=4,outline=TRUE,fsize=c(0.7,0.9),legend=50)

## -- Or we can show the full posterior sample of changes of each type on the
## -- tree. These are not changes from any single stochastic character history,
## -- but changes in our full posterior sample of such histories.

dotTree(eel.tree,fmode,colors=setNames(c("red","blue"),c("bite","suction")),
    fsize=0.7,ftype="i",legend=FALSE)
add.simmap.legend(x=0,y=-4,colors=setNames(c("red","blue"),c("bite",
    "suction")),prompt=FALSE,vertical=FALSE,shape="circle")
nulo<-sapply(eel.trees,markChanges,sapply(setNames(c("red","blue"),
    c("bite","suction")),make.transparent,0.05))
add.simmap.legend(colors=sapply(setNames(setNames(c("red","blue"),
    c("bite","suction"))[2:1],
    c("bite->suction","suction->bite")),
    make.transparent,0.1),prompt=FALSE,x=50,y=-4,vertical=FALSE)

## -- In addition to these methods for either continuous or discretely 
## -- measured character traits, we can easily combine the methods.

## -- For instance, in the following I will overlay the posterior density
## -- from stochastic mapping on a dotTree plot of the continuous trait of
## -- overall body size.

bsize<-setNames(eel.data[,2],rownames(eel.data))
dotTree(eel.tree,bsize,length=10,fsize=0.7,lwd=7,ftype="i")
text(x=28,y=-1.8,"body size (cm)",pos=4)
plot(obj$tree,colors=obj$cols,add=TRUE,ftype="off",lwd=5,
    xlim=get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim,
    ylim=get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim)
add.color.bar(leg=50,cols=obj$cols,title="PP(state=bite)",
    prompt=FALSE,x=70,y=-2)

## -- Alternatively, we could do this as two facing plots - one of densityMap
## -- and the other of contMap styles:

layout(matrix(1:3,1,3),widths=c(0.39,0.22,0.39))
plot(obj,lwd=6,ftype="off",legend=60,outline=TRUE,fsize=c(0,1.2))
plot.new()
plot.window(xlim=c(-0.1,0.1),
    ylim=get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim)
par(cex=0.6)
text(rep(0,length(obj$tree$tip.label)),1:Ntip(obj$tree),
    gsub("_"," ",obj$tree$tip.label),font=3)
plot(setMap(contMap(eel.tree,bsize,plot=FALSE),invert=TRUE),lwd=6,outline=TRUE,
    direction="leftwards",ftype="off",legend=60,fsize=c(0,1.2),
    leg.txt="body size (cm)")

## -- As another example here is a stochastically mapped discrete character
## -- history overlain on a traitgram plot of body size. At nodes I will add
## -- the posterior probabilties from stochastic mapping.

phenogram(eel.trees[[1]],bsize,lwd=3,colors=
    setNames(c("blue","red"),c("suction","bite")),
    spread.labels=TRUE,spread.cost=c(1,0),fsize=0.6,
    ftype="i")

## Optimizing the positions of the tip labels...

add.simmap.legend(colors=setNames(c("blue","red"),c("suction","bite")),
    prompt=FALSE,shape="circle",x=0,y=250)
obj<-summary(eel.trees)
nodelabels(pie=obj$ace,piecol=setNames(c("red","blue"),colnames(obj$ace)),
    cex=0.6)
tiplabels(pie=to.matrix(fmode[eel.tree$tip.label],colnames(obj$ace)),
    piecol=setNames(c("red","blue"),colnames(obj$ace)),cex=0.4)

## -- These are just a few examples. There are many other ways in which 
## -- discrete & continuous character methods can be combined

## -- I'm now going to demonstrate a handful of other phytools plotting methods.

## -- Firstly, phytools can be used to project a tree onto a geographic map.
## -- This can be done using linking lines from the tips of a square phylogram
## -- or directly.

lat.long<-read.csv("latlong.csv",row.names=1)
obj<-phylo.to.map(tree,lat.long,plot=FALSE)

## objective: 110
## objective: 110
## objective: 110
## objective: 110
## objective: 110
## objective: 110
## objective: 110
## objective: 110

## objective: 100
## objective: 100

## objective: 98
## objective: 98
## objective: 98
## objective: 98
## objective: 98
## objective: 98

## objective: 96

## objective: 92
## objective: 92

## objective: 86
## objective: 86
## objective: 86
## objective: 86
## objective: 86
## objective: 86

plot(obj,type="phylogram",asp=1.3,mar=c(0.1,0.5,3.1,0.1))

plot(obj,type="direct",asp=1.3,mar=c(0.1,0.5,3.1,0.1))

## -- phytools can be used to create a co-phylogenetic plot in which two
## -- trees are graphed in a facing fashion, with linking lines between
## -- corresponding tips. The following uses data from Lopez-Vaamonde et al. 
## -- (2001; MPE).

trees<-read.tree("cophylo.tre")
Pleistodontes<-trees[[1]]
Sycoscapter<-trees[[2]]
assoc<-read.csv("assoc.csv")
obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
    print=TRUE)

## Rotating nodes to optimize matching...

## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778

## objective: 357.777777777778
## objective: 357.777777777778

## objective: 345.111111111111
## objective: 345.111111111111
## objective: 345.111111111111

## objective: 322.777777777778

## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856

## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778

## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856

## Done.

plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8)

## -- This method is fairly flexible. For instance, we can use add tip, node,
## -- or edge labels to each tree, or we can us different colors for different 
## -- edges.

plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col="grey",fsize=0.8)
nodelabels.cophylo(which="left",frame="circle",cex=0.8)
nodelabels.cophylo(which="right",frame="circle",cex=0.8)

left<-rep("black",nrow(obj$trees[[1]]$edge))
nodes<-getDescendants(obj$trees[[1]],30)
left[sapply(nodes,function(x,y) which(y==x),y=obj$trees[[1]]$edge[,2])]<-
    "blue"
right<-rep("black",nrow(obj$trees[[2]]$edge))
nodes<-getDescendants(obj$trees[[2]],25)
right[sapply(nodes,function(x,y) which(y==x),y=obj$trees[[2]]$edge[,2])]<-
    "blue"
edge.col<-list(left=left,right=right)
plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8,
    edge.col=edge.col)

## -- phytools contains a method for plotting a posterior sample of trees from
## -- Bayesian analysis. This is more or less equivalent to the stand-alone
## -- program densiTree. Here's a demo:

trees<-read.tree("density.tre")
densityTree(trees,type="cladogram",nodes="intermediate")

## -- We can also just represent topological incongruence among trees:

densityTree(trees,use.edge.length=FALSE,type="phylogram",nodes="inner")

## -- If we want to just compare two alternative time trees for the same set
## -- of taxa, we could try compare.chronograms. Here is a demo with simulated
## -- trees.

tree<-rtree(n=26,tip.label=LETTERS)
t1<-force.ultrametric(tree,"nnls")
tree$edge.length<-runif(n=nrow(tree$edge))
t2<-force.ultrametric(tree,"nnls")
compare.chronograms(t1,t2)

## -- phytools has an interactive node labeling method. Here's a demo using a
## -- phylogeny of vertebrates (in non-interactive mode)

vertebrates<-read.tree("vertebrates.tre")
plotTree(vertebrates,xlim=c(-24,450))
labels<-c("Cartilaginous fish",
    "Ray-finned fish",
    "Lobe-finned fish",
    "Anurans",
    "Reptiles (& birds)",
    "Birds",
    "Mammals",
    "Eutherians")
nodes<-labelnodes(text=labels,node=c(24,27,29,32,34,35,37,40),
    shape="ellipse",cex=0.8,interactive=FALSE)

## -- Another method in phytools allows us to easily label only a subset of tips
## -- in the tree. This can be handy for very large trees in which we are only
## -- interested in labeling a handful of taxa.

plot(anoletree,cols,xlim=c(0,9),ylim=c(-4,Ntip(anoletree)),ftype="off")
add.simmap.legend(colors=cols,prompt=FALSE,x=0,y=-4,vertical=FALSE)
text<-sample(anoletree$tip.label,20)
tips<-sapply(text,function(x,y) which(y==x),y=anoletree$tip.label)
linklabels(text,tips,link.type="curved")

## -- Here is a recent method for adding circular clade labels to a 
## -- plotted fan style tree.

obj<-setMap(contMap(anole.tree,svl,plot=FALSE),invert=TRUE)
plot(obj,type="fan",ftype="off",lwd=c(2,6),outline=FALSE,xlim=c(-1.3,1.3),
    legend=0.8,leg.txt="log(SVL)")
nodes<-c(105,117,122,135,165,169,68,172,176,183,184)
labels<-paste("Clade",LETTERS[1:length(nodes)])
for(i in 1:length(nodes))
    arc.cladelabels(tree=obj$tree,labels[i],nodes[i],
        orientation=if(labels[i]%in%c("Clade G","Clade J"))
        "horizontal" else "curved",mark.node=FALSE)

## -- Finally, the following is a method to add a geological timescale to a plotted
## -- tree.

tree<-pbtree(b=0.03,d=0.01,n=200)
h<-max(nodeHeights(tree))
plotTree(tree,plot=FALSE)
obj<-geo.legend(alpha=0.3,cex=1.2,plot=FALSE)
obj$leg<-h-obj$leg
plotTree(tree,ftype="off",ylim=c(-0.2*Ntip(tree),Ntip(tree)),lwd=1)
geo.legend(leg=obj$leg,colors=obj$colors,cex=1.2)

That's it!

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!