Online practical materials for plotting methods book chapter

I have been working on my long overdue 'online practical materials' for the Garamszegi book chapter. My intention is to post code so that readers could reproduce all the figures of the chapter. My chapter is available here. My OPM draft is below. The final OPM will be posted here.

Figures for “Chapter 4: Graphical Methods for Visualizing Comparative Data on Phylogenies”

In the following twelve sections I will describe and provide code that can be used to reproduce the twelve figures of my chapter in Modern Phylogenetic Comparative Methods and their Application in Evolutionary Biology: Concepts and Practice, L´szló Zsolt Garamszegi Ed. Please contact me if you have any questions about this material.

Figure 4.1:

Figure 4.1 is just a simple illustration of two different tree plotting styles that are employed repeatedly throughout the chapter: a right-facing square phylogram; and a circular (fan) phylogram. This exercise requires the file Centrarchidae.tre. Here's how the plots were produced:

library(phytools)
packageVersion("phytools")

## [1] '0.4.38'

text<-"((L:2.29,(K:1.33,(J:1.03,I:1.03):0.3):0.96):0.19,(((H:1.05,(G:0.75,(F:0.73,E:0.73):0.01):0.3):0.71,(D:1.07,(C:0.96,B:0.96):0.11):0.69):0.19,A:1.95):0.53);"
tree<-read.tree(text=text)
rm(text)
centrarchidae<-read.tree(file="Centrarchidae.tre")
par(mfcol=c(1,2))
plotTree(tree,mar=rep(1,4),fsize=1.5)
title(main="(a)",adj=0,cex.main=2)
plotTree(centrarchidae,type="fan",mar=rep(1,4))

## setEnv=TRUE for this type is experimental. please be patient with bugs

title(main="(b)",adj=0,cex.main=2)

Figure 4.2:

Figure 4.2 illustrates, step-by-step, the algorithm for drawing a rightward-facing, square phylogram. Here is the code illustrating how this figure was produced:

## preliminaries
n<-length(tree$tip.label)
# reorder cladewise to assign tip positions
cw<-reorder(tree,"cladewise")
y<-vector(length=n+cw$Nnode)
y[cw$edge[cw$edge[,2]<=n,2]]<-1:n
# reorder pruningwise for post-order traversal
pw<-reorder(tree,"pruningwise")
nn<-unique(pw$edge[,1])
# compute vertical position of each edge
for(i in 1:length(nn)){
    yy<-y[pw$edge[which(pw$edge[,1]==nn[i]),2]]
    y[nn[i]]<-mean(range(yy))
}
# compute start & end points of each edge
X<-nodeHeights(cw)

## end preliminaries

matrix(c(1,2,3,6,5,4),3,2,byrow=FALSE)->ii
matrix(1:6,3,2,byrow=FALSE)->ii
matrix(1:6,3,2,byrow=TRUE)->ii

layout(ii)

## plot (a)

# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(a)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")

## NULL

for(i in 1:n)
    text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)

## plot (b)

# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(b)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")

## NULL

for(i in 1:n)
    text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)

for(i in 2:tree$Nnode+length(tree$tip.label)){
    lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),col="gray",
        lty="dashed")
    text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
        paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,i)<=length(cw$tip)]],
        collapse=","),col="gray",pos=4,cex=2)
}

## plot (c)

# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(c)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",
    lty="dashed")

## NULL

for(i in 1:n)
    text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)

for(i in 2:tree$Nnode+length(tree$tip.label)){
    lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),col="gray",
        lty="dashed")
    text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
        paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,i)<=length(cw$tip)]],
        collapse=","),col="gray",pos=4,cex=2)
}

for(i in 1:nrow(X))
    points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)

## plot (d)

# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(d)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")

## NULL

for(i in 1:n)
    text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)

for(i in 2:tree$Nnode+length(tree$tip.label)){
    lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),
        col="gray",lty="dashed")
    text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
        paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,
        i)<=length(cw$tip)]],collapse=","),col="gray",pos=4,cex=2)
}

# plot horizontal edges
for(i in 1:nrow(X))
    lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=2,lend=2)

for(i in 1:nrow(X))
    points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)

## plot (e)

# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(e)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")

## NULL

## abline(h=1:max(y),col="gray",lty="dashed")

for(i in 1:n)
    text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray")

for(i in 2:tree$Nnode+length(tree$tip.label)){
    lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),
        col="gray",lty="dashed")
    text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
        paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,
        i)<=length(cw$tip)]],collapse=","),col="gray",pos=4,cex=2)
}

# plot horizontal edges
for(i in 1:nrow(X))
    lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=2,lend=2)

# plot vertical relationships
for(i in 1:tree$Nnode+n)
    lines(X[which(cw$edge[,1]==i),1],range(y[cw$edge[which(cw$edge[,1]==i),
        2]]),lwd=2,lend=2)

for(i in 1:nrow(X))
    points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)

## plot (f)

# open & size a new plot
plot.new(); par(mar=c(1.1,1.1,1.5,0.1))
title(main="(f)",adj=0,cex.main=3)
plot.window(xlim=c(0,1.05*max(X)),ylim=c(1,max(y)+1))
axis(1,at=c(0,max(X)),labels=FALSE); axis(2,at=1:max(y),labels=FALSE)
apply(cbind(1:max(y),1:max(y)),1,lines,x=c(0,max(X)),col="gray",lty="dashed")

## NULL

## abline(h=1:max(y),col="gray",lty="dashed")

for(i in 1:n)
    text(max(X),y[i],cw$tip.label[i],pos=4,offset=0.3,col="gray",cex=2)

for(i in 2:tree$Nnode+length(tree$tip.label)){
    lines(x=c(0,max(X)),y=rep(y[cw$edge[which(tree$edge[,2]==i),2]],2),
        col="gray",lty="dashed")
    text(X[which(cw$edge[,1]==i),1],y[cw$edge[which(tree$edge[,2]==i),2]],
        paste(cw$tip.label[getDescendants(cw,i)[getDescendants(cw,
        i)<=length(cw$tip)]],collapse=","),col="gray",pos=4,cex=2)
}

# plot horizontal edges
for(i in 1:nrow(X))
    lines(X[i,],rep(y[cw$edge[i,2]],2),lwd=2,lend=2)

# plot vertical relationships
for(i in 1:tree$Nnode+n)
    lines(X[which(cw$edge[,1]==i),1],range(y[cw$edge[which(cw$edge[,1]==i),
        2]]),lwd=2,lend=2)

for(i in 1:nrow(X))
    points(X[i,],rep(y[cw$edge[i,2]],2),pch=21,bg="gray",cex=2)


# plot tip labels
for(i in 1:n)
    text(X[which(cw$edge[,2]==i),2],y[i],tree$tip.label[i],pos=4,offset=0.3,
        font=2,cex=2)

Figure 4.3:

Figure 4.3 shows a single stochastic character map with the aggregate posterior probabilities at all nodes from 100 stochastic character maps overlain. Here is how it was created:

# load anole 'ecomorph' tree
data(anoletree)
# this is a stochastic mapped tree - pull tip states
x<-getStates(anoletree,"tips")
# conduct stochastic mapping
trees<-make.simmap(anoletree,x,model="ER",nsim=100,message=FALSE)
# get the states at ancestral nodes
aa<-describe.simmap(trees,plot=FALSE)
# set colors to be used for mapping
# cols<-setNames(palette()[1:6],sort(unique(x)))
cols<-setNames(c("blue","gold1","chartreuse3","saddlebrown",
    "orangered","purple"),sort(unique(x)))
# plotSimmap
plotSimmap(anoletree,cols,type="fan",lwd=3,fsize=0.8,ftype="i")

## setEnv=TRUE for this type is experimental. please be patient with bugs

# add node labels
par(fg="transparent")
nodelabels(pie=aa$ace,piecol=cols,cex=0.4)
# add legend
par(fg="black")
add.simmap.legend(colors=cols,fsize=0.8,x=-max(nodeHeights(anoletree)),
    y=-max(nodeHeights(anoletree)),prompt=FALSE)

Figure 4.4:

Figure 4.4 has two parts. The first part shows a posterior density map of feeding mode on the tree of elopomorph eels, and the second part shows the variance among stochastic character maps for the same tree & trait. This exercise requires the files elopomorph.tre and elophomorph.csv. Here is how that figure was created:

tree<-read.tree("elopomorph.tre")
x<-as.matrix(read.csv("elopomorph.csv",row.names=1))[,1]
mtrees<-make.simmap(tree,x,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##                 biting suction feeding
## biting          -19.02           19.02
## suction feeding  19.02          -19.02

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##          biting suction feeding 
##             0.5             0.5

## Done.

dMap<-densityMap(mtrees,states=c("suction feeding","biting"),plot=FALSE)

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

## create 'variance' map
vMap<-dMap
aa<-colnames(vMap$tree$mapped.edge)
aa<-colnames(vMap$tree$mapped.edge)
tt<-vMap$tree
# collapse each prob with the same variance into each other
for(i in 0:499){
    oo<-c(as.character(i),as.character(1000-i))
    nn<-as.character(i)
    chk<-oo%in%aa
    if(sum(chk)==2) tt<-mergeMappedStates(tt,oo,nn)
    else if(chk[2])
    for(j in 1:length(tt$maps)) 
        names(tt$maps[[j]])[which(names(tt$maps[[j]])==oo[2])]<-oo[1]
}
# compute the variances
p<-as.numeric(names(vMap$cols))/(length(vMap$cols)-1)
v<-(p*(1-p))[1:501]
vMap$tree<-tt
# this is a hack
vMap$cols<-setNames(as.vector(rbind(grey(v/0.25), grey(v/0.25))),
    as.vector(rbind(0:500,0:500)))
vMap$lims<-c(0,0.25)
class(vMap)<-"contMap"

par(mfrow=c(1,2))

plot(dMap,fsize=c(0.7,0.9),outline=TRUE,legend=0.05,mar=c(0.1,0.5,2.1,0.1))
title(main="(a)",adj=0)

par(col="white")
plotTree(tree,fsize=0.7,lwd=3+2,offset=0.2+0.2/3,ftype="i",ylim=c(-6.32,62),
    mar = c(0.1,0.5,2.1,0.1))
par(col="black")
plotSimmap(vMap$tree,vMap$cols,lwd=3,fsize=0.7,ftype="i",add=TRUE,
    ylim=c(-6.32,62),mar=c(0.1,0.5,2.1,0.1))

add.color.bar(0.05,vMap$cols,title="variance",lims=as.numeric(vMap$lims),
    digits=2,x=0,y=1-0.08*(62-1),lwd=3,fsize=0.9,prompt=FALSE)
title(main="(b)",adj=0)

Figure 4.5:

Figure 4.5 is a three-part illustration of one & two-dimensional traitgrams with simulated data. Since I didn't save the seed, results will vary. Here is how it was created:

tree<-pbtree(n=26)
tree$tip.label<-LETTERS
X<-sim.corrs(tree,matrix(c(1,0.7,0.7,1),2,2))
par(mfcol=c(1,3))
par(mar=c(4.1,4.1,4.1,2.1))
phenogram(tree,X[,1],spread.labels=TRUE)
title(main="(a)",adj=0,cex=3)
fancyTree(tree,"traitgram3d",X=X,method="static",angle=-30)
title(main="(b)",adj=0,cex=3)
par(mar=c(4.1,4.1,4.1,2.1))
fancyTree(tree,"phenogram95",x=X[,1],link=0.1*max(nodeHeights(tree)))
title(main="(c)",adj=0,cex=3)

Figure 4.6:

Figure 4.6 is a projection of the phylogeny into a two dimensional morphospace (aka. a phylomorphospace) with the mapping of a multistate discrete ecological trait overlain. The tree and data are for Caribbean Anolis lizards. This exercise requires the file anole.data.csv. Here is how it was created:

X<-read.csv("anole.data.csv",header=T,row.names=1)
X<-as.matrix(X)
rownames(X)<-paste("Anolis_",rownames(X),sep="")
data(anoletree)
X<-X[anoletree$tip.label,]
X<-cbind(X[,1],phyl.resid(anoletree,X[,1],X[,2:20])$resid,phyl.resid(anoletree,
    X[,21],X[,22])$resid)
PC<-phyl.pca(anoletree,X)$S
PC<-cbind(-PC[,1],PC[,2])
colnames(PC)<-c("PC1 (limb length)","PC2 (size)")
par(mar=c(5.1,4.1,2.1,2.1))
colors<-setNames(c("blue","gold1","chartreuse3","saddlebrown","orangered","purple"),
    sort(unique(getStates(anoletree,"tips"))))
phylomorphospace(anoletree,PC[,c(1,2)],label="off",node.by.map=TRUE,colors=colors)
add.simmap.legend(colors=colors,prompt=FALSE,x=-1.5,y=1.1)

Figure 4.7:

Figure 4.7 is a stochastic phylogeny and projection of that phylogeny into bivariate morphospace with a continuous color gradient showing time since the root overlain. The purpose of this visualization is to also convey temporal information in a phylomorphospace projection. This plot was created as follows (although I neglected to save the seed in the manuscript version, so this results will differ from the published one):

# simulate a tree
set.seed(1)
tree<-pbtree(n=26,scale=100)
tree$tip.label<-LETTERS[26:1]
# simulate data for our phylomorphospace
X<-fastBM(tree,nsim=2)
colnames(X)<-paste("V",1:2,sep="")
# now add the era map
tree<-make.era.map(tree,0:99)
# choose the color scale
colors<-rainbow(100,start=0,end=0.7)
colors<-colorRampPalette(c("blue", "red"))(100)
colors<-colorRampPalette(c("red", "orange", "blue"),space="rgb")(100)
names(colors)<-1:100
par(mfcol=c(1,2))
# check it
plotSimmap(tree,colors,pts=FALSE,lwd=6,mar=c(5.1,2.1,4.1,1.1),ftype="i")
title(main="(a)",adj=0)
# plot phylomorphospace
par(mar=c(5.1,4.1,4.1,1.1))
phylomorphospace(tree,X,colors=colors,node.by.map=TRUE,lwd=4,node.size=c(0.8,1.3),
    label="horizontal")
add.color.bar(cols=colors,leg=10,title="time since the root",subtitle="",
    lims=c(0,max(nodeHeights(tree))),prompt=FALSE,x=-7.5,y=16)
title(main="(b)",adj=0)

Figure 4.8:

Figure 4.8 shows observed and reconstructed body size (snout-to-vent length, SVL) mapped onto a phylogeny of Greater Antillean Anolis lizard species. This exercise also requires the file anole.data.csv.

X<-read.csv("anole.data.csv",header=TRUE,row.names=1)
data(anoletree)
svl<-as.matrix(X)[,1]
names(svl)<-paste("Anolis_",names(svl),sep="")
svl<-svl[anoletree$tip.label]
obj<-contMap(anoletree,svl,plot=FALSE)
plot(obj,type="fan",legend=5)

Figure 4.9:

Figure 4.9 shows a simlated four trait “phylogenetic scatterplot matrix” continuous character color mappings of each trait on the diagonal, and bivariate phylomorphospaces in off-diagonal positions. Once again, this result uses simulation so results will vary.

tree<-pbtree(n=26,tip.label=LETTERS[26:1])
V<-matrix(c(1,0.5,0.2,0.0,
        0.5,1,0.7,0.2,
        0.2,0.7,1,0.1,
        0.0,0.2,0.1,1),4,4)
V

##      [,1] [,2] [,3] [,4]
## [1,]  1.0  0.5  0.2  0.0
## [2,]  0.5  1.0  0.7  0.2
## [3,]  0.2  0.7  1.0  0.1
## [4,]  0.0  0.2  0.1  1.0

X<-sim.corrs(tree,V)
fancyTree(tree,type="scattergram",X=X)

Figure 4.10:

Figure 4.10 shows in panel (a) a traitgram with a discrete character history overlain. Data were simulated with a high rate of evolution on the red branches. Panel (b) shows a posterior density map from stochastic character mapping projected onto a phylomorphospace plot. In this case, data were simulated with a higher evolutionary correlation on blue than red branches.

set.seed(4)
t1<-pbtree(n=26,scale=1)
t1$tip.label<-LETTERS
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-0:1
t1<-sim.history(t1,Q,anc="0")

## Done simulation(s).

colors<-setNames(c("blue","red"),0:1)

par(lend=2)
x<-sim.rates(t1,c(1,10))

## names absent from sig2: assuming same order as $mapped.edge

par(mfcol=c(1,2))
phenogram(t1,x,colors=colors,spread.labels=TRUE,lwd=3,spread.cost=c(1,0),ftype="i")
title(main="(a)",adj=0)

V<-setNames(
    list(matrix(c(1,0.8,0.8,1),2,2),matrix(c(2,0,0,2),2,2)),
    0:1)
X<-sim.corrs(t1,V)
colnames(X)<-paste("V",1:2,sep="")
trees<-make.simmap(t1,t1$states,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##        0      1
## 0 -2.683  2.683
## 1  2.683 -2.683

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##   0   1 
## 0.5 0.5

## Done.

xx<-densityMap(trees,plot=FALSE)

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

phylomorphospace(xx$tree,X,colors=xx$cols,lwd=3,node.by.map=TRUE,label="horizontal")
title(main="(b)",adj=0)

Figure 4.11:

Figure 4.11 shows two different projections of a tree onto a geographic map using simulated data. The following shows how this was created:

# simulate tree & data
tree<-pbtree(n=26,scale=100)
tree$tip.label<-LETTERS[26:1]
lat<-fastBM(tree,sig2=10,bounds=c(-90,90))
long<-fastBM(tree,sig2=80,bounds=c(-180,180))
# now plot projection
xx<-phylo.to.map(tree,cbind(lat,long),plot=FALSE)

## objective: 216
## objective: 206
## objective: 206
## objective: 206
## objective: 206
## objective: 206
## objective: 202
## objective: 200
## objective: 200
## objective: 198
## objective: 198
## objective: 172
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154
## objective: 154

layout(matrix(c(1,2),2,1),heights=c(0.61,0.39))
plot(xx,type="phylogram",asp=1.3,mar=c(0.1,0.5,3.1,0.1))
title(main="(a)",adj=0)
plot(xx,type="direct",asp=1.3,mar=c(0.1,0.5,3.1,0.1))
title(main="(b)",adj=0)

Figure 4.12:

Figure 4.12 helps illustrate the structure of an object of class "phylo" in memory in R. Here is how it was created:

tree<-list()
class(tree)<-"phylo"
tree$edge<-matrix(c(6,7,
            7,1,
            7,8,
            8,9,
            9,2,
            9,3,
            8,4,
            6,5),8,2,byrow=TRUE)
tree$Nnode<-4
tree$tip.label<-LETTERS[1:5]
tree$edge.length<-c(4,8,5,1,2,2,3,12)

library(plotrix)
plotTree(tree,setEnv=TRUE,offset=1,fsize=1.5,ftype="i",ylim=c(0.5,5.5))
lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
boxed.labels(x=lastPP$xx,y=lastPP$yy,1:9,bg=grey(level=0.45),cex=1.6,xpad=2.0,
    ypad=1.5)

Liam J. Revell, Nov. 26 2014.

Pruning trees to one member per genus, and to one descendant for each clade younger than a particular age

The following request was recently posted to the R-sig-phylo email list:

“(1) I'd like to drop all but one of the tips in each genus (i.e. to generate a genus-level phylogeny). Given that all genera are monophyletic, it shouldn't matter which species I pick.

(2) I'd like to establish an era (e.g. 25 Mya) and leave only one representative of the clades there were present at that point in time.”

I submitted my replies (1, 2) but I'll repeat them here anyway.

First, let's simulate scenario (1) in which we have a species-level tree with one or more species per genus, and want to reduce the tree to a genus tree:

library(phytools)
## here are our species names
tips<-c("Abc_def","Abc_ghi","Def_ghi","Def_jkl",
    "Ghi_jkl","Ghi_mno","Jkl_mno","Jkl_pqr")
tree<-stree(n=8,type="balanced",tip.label=tips)
plotTree(tree,ftype="i")

Next, we'll get a list of all genus names, and prune the tree of all but one member of each genus:

## get a list of all genera
tips<-tree$tip.label
genera<-unique(sapply(strsplit(tips,"_"),function(x) x[1]))
## here are our genera
genera

## [1] "Abc" "Def" "Ghi" "Jkl"

## now drop all but one of each
ii<-sapply(genera,function(x,y) grep(x,y)[1],y=tips)
tree<-drop.tip(tree,setdiff(tree$tip.label,tips[ii]))
plotTree(tree,ftype="i")

Finally, we can rename our tips to be the genus names only:

tree$tip.label<-sapply(strsplit(tree$tip.label,"_"),function(x) x[1])
plotTree(tree,ftype="i")

That's it.

Now, let's move on to part (2). For this, I'll simulate a true tree of arbitrary total depth 100.

tree<-pbtree(n=100,scale=100)
plotTree(tree,fsize=0.6)

Now, let's prune each clade younger than 25 units of time old to include only 1 descendant:

## get all node heights
H<-nodeHeights(tree)
## time from the root
t<-max(H)-25
## identify all edges that cross 25 mybp
h1<-which(H[,1]<t)
h2<-which(H[,2]>t)
ii<-intersect(h1,h2)
## all daughter nodes of those edges
nodes<-tree$edge[ii,2]
## internal phytools function
getDescendants<-phytools:::getDescendants
## find all descendants from each edge
tips<-lapply(nodes,getDescendants,tree=tree)
## find all tips to keep
tips<-tree$tip.label[sapply(tips,function(x,y) x[x<=Ntip(y)][1],y=tree)]
## drop all the others
tree<-drop.tip(tree,setdiff(tree$tip.label,tips))
plotTree(tree)
lines(c(t,t),par()$usr[3:4],lty="dashed",col="red")

Cool. It worked!

That's all for now.

Alternative implementations of fitContinuous

A colleague contacted me today about “alternative implementations” of geiger's fitContinuous function that he could use to check parameter estimation, and maybe to speed calculations as fitContinuous can be somewhat slow (although is evidently quite robust now) for some models.

Well, I don't know about speeding things up - but it is pretty straightforward to write a wrapper function that combines phytools brownie.lite single-rate model with geiger rescale.phylo to fit many of the models of fitContinuous. The following is a simple demo:

## load packages
library(phytools)
library(geiger)

## helper function to get the root node number
getRoot<-function(tree) Ntip(tree)+1

## here is our fitContinuous lite function
fitCont<-function(tree,x,model="BM",interval=NULL){
    if(model!="BM"){
        lk<-function(par,tree,x,model) brownie.lite(paintSubTree(rescale(tree,model,par),getRoot(tree),"1"),x)$logL1
        oFit<-optimize(lk,interval,tree=tree,x=x,model=model,maximum=TRUE)
        cFit<-brownie.lite(paintSubTree(rescale(tree,model,oFit$maximum),getRoot(tree),"1"),x)
        obj<-list(par=oFit$maximum,model=model,sig2=cFit$sig2.single,a=cFit$a.single,logLik=oFit$objective)
    } else {
        fit<-brownie.lite(paintSubTree(tree,getRoot(tree),"1"),x)
        obj<-list(model=model,sig2=fit$sig2.single,a=fit$a.single,logLik=fit$logL1)
    }
    obj
}

OK, now let's try it for simulated data:

## simulate tree
tree<-pbtree(n=100,scale=1)
## simulate data under model="BM"
x<-fastBM(tree)
fitContinuous(tree,x)

## GEIGER-fitted comparative model of continuous data
##  fitted 'BM' model parameters:
##  sigsq = 0.940924
##  z0 = -0.186971
## 
##  model summary:
##  log-likelihood = -58.714396
##  AIC = 121.428793
##  AICc = 121.552504
##  free parameters = 2
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 1.00
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

fitCont(tree,x)

## $model
## [1] "BM"
## 
## $sig2
## [1] 0.9409
## 
## $a
## [1] -0.187
## 
## $logLik
## [1] -58.71

## simulate data under model="lambda"
x<-fastBM(rescale(tree,model="lambda",0.5))
fitContinuous(tree,x,model="lambda")

## GEIGER-fitted comparative model of continuous data
##  fitted 'lambda' model parameters:
##  lambda = 0.729756
##  sigsq = 1.217453
##  z0 = 0.254493
## 
##  model summary:
##  log-likelihood = -122.552971
##  AIC = 251.105941
##  AICc = 251.355941
##  free parameters = 3
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 0.32
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

fitCont(tree,x,model="lambda",interval=c(0,1))

## $par
## [1] 0.7298
## 
## $model
## [1] "lambda"
## 
## $sig2
## [1] 1.217
## 
## $a
## [1] 0.2545
## 
## $logLik
## [1] -122.6

## simulate data under model="OU"
x<-fastBM(rescale(tree,model="OU",0.4))
fitContinuous(tree,x,model="OU")

## GEIGER-fitted comparative model of continuous data
##  fitted 'OU' model parameters:
##  alpha = 0.536268
##  sigsq = 1.171551
##  z0 = 0.216378
## 
##  model summary:
##  log-likelihood = -61.980067
##  AIC = 129.960133
##  AICc = 130.210133
##  free parameters = 3
## 
## Convergence diagnostics:
##  optimization iterations = 100
##  failed iterations = 0
##  frequency of best fit = 0.60
## 
##  object summary:
##  'lik' -- likelihood function
##  'bnd' -- bounds for likelihood search
##  'res' -- optimization iteration summary
##  'opt' -- maximum likelihood parameter estimates

fitCont(tree,x,model="OU",interval=c(0,1))

## $par
## [1] 0.5363
## 
## $model
## [1] "OU"
## 
## $sig2
## [1] 1.172
## 
## $a
## [1] 0.2164
## 
## $logLik
## [1] -61.98

Running system.time with any of these examples (except for model="BM") does not show any improvement of this ostensibly lighter implementation - but this is probably due to lots of unnecessary internals in brownie.lite for fitting a multi-rate model (which is not actually done here). Improving on this is thus pretty simple (but I will not do so here, nor now).

That's it.

Optimizing tree transformations for multiple discrete characters

A recent R-sig-phylo subscriber asked:

“I have a tree and discrete data (number of gene copies, for many genes) and would like to use the fitDiscrete function in geiger, or something similar. However, I would like to estimate the parameters given all of the datasets, not just with the data for each gene. For instance, if I was using the "delta” model to vary rates across the tree, I would like this delta value to reflect some sort of summary value across all datasets. Does anyone have an idea as to how this could be accomplished or perhaps point me in the right direction?“

Brian O'Meara pointed out that it would be straighforward to write a likelihood function that just summed the log-likelihoods across characters for different tree transformations and then wrapped the function in an optimizer. The following illustrates this explicitly.

First load packages - for this we need ape, phytools (just for simulation here), and geiger (for our tree rescaling):

## load packages
library(ape)
library(phytools)

## Loading required package: maps

library(geiger)

Now let's simulate some data that would be analogous to what we'd load from file in the empirical case:

## simulate tree & data
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
tree<-pbtree(n=100,scale=1)
X<-replicate(10,sim.history(rescale(tree,model="lambda",lambda=0.5),Q)$states)

## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).
## Done simulation(s).

## this is what our data look like:
X

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## t64  "a"  "a"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  "a"  
## t97  "a"  "b"  "b"  "a"  "a"  "b"  "b"  "b"  "a"  "a"  
## t98  "b"  "a"  "a"  "b"  "b"  "b"  "b"  "b"  "b"  "a"  
## t43  "b"  "b"  "b"  "a"  "b"  "a"  "a"  "b"  "a"  "b"  
## t44  "b"  "a"  "a"  "b"  "b"  "b"  "a"  "a"  "a"  "b"  
## t23  "b"  "a"  "b"  "a"  "b"  "a"  "a"  "a"  "b"  "a"  
## t85  "b"  "a"  "b"  "b"  "a"  "a"  "b"  "b"  "b"  "b"  
## t86  "b"  "a"  "b"  "b"  "b"  "b"  "a"  "b"  "a"  "a"  
## t68  "a"  "a"  "b"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  
## t83  "b"  "a"  "b"  "b"  "b"  "a"  "a"  "b"  "a"  "b"  
## t84  "a"  "a"  "b"  "a"  "b"  "a"  "b"  "a"  "b"  "b"  
## t2   "a"  "a"  "a"  "a"  "b"  "b"  "a"  "b"  "b"  "a"  
## t3   "b"  "a"  "a"  "a"  "b"  "a"  "a"  "a"  "a"  "a"  
## t91  "b"  "b"  "b"  "a"  "a"  "a"  "b"  "a"  "a"  "a"  
## t92  "a"  "a"  "a"  "a"  "a"  "a"  "b"  "a"  "a"  "a"  
## t45  "b"  "a"  "b"  "a"  "a"  "b"  "a"  "a"  "a"  "b"  
## t46  "b"  "a"  "b"  "b"  "b"  "a"  "a"  "b"  "a"  "a"  
## t58  "b"  "b"  "a"  "a"  "a"  "a"  "a"  "a"  "a"  "a"  
## t59  "a"  "a"  "b"  "b"  "b"  "b"  "a"  "b"  "b"  "b"  
## t21  "b"  "b"  "b"  "b"  "b"  "a"  "a"  "a"  "b"  "b"  
## t60  "b"  "a"  "a"  "a"  "b"  "b"  "b"  "a"  "b"  "b"  
## t89  "a"  "a"  "b"  "a"  "b"  "b"  "a"  "a"  "a"  "b"  
## t90  "b"  "b"  "a"  "b"  "a"  "a"  "a"  "a"  "b"  "b"  
## t36  "b"  "a"  "b"  "a"  "b"  "b"  "b"  "b"  "b"  "a"  
## t27  "b"  "b"  "b"  "b"  "b"  "a"  "b"  "a"  "a"  "b"  
## t28  "a"  "b"  "b"  "b"  "b"  "b"  "a"  "a"  "b"  "b"  
## t29  "b"  "b"  "b"  "a"  "a"  "b"  "a"  "a"  "b"  "b"  
## t37  "b"  "b"  "b"  "b"  "b"  "a"  "a"  "a"  "b"  "b"  
## t38  "b"  "a"  "a"  "a"  "a"  "b"  "a"  "a"  "a"  "b"  
## t16  "a"  "a"  "b"  "b"  "a"  "a"  "a"  "b"  "a"  "a"  
## t17  "a"  "a"  "a"  "b"  "a"  "b"  "b"  "a"  "a"  "b"  
## t15  "b"  "a"  "a"  "a"  "a"  "b"  "b"  "a"  "b"  "b"  
## t47  "b"  "a"  "b"  "a"  "b"  "b"  "b"  "a"  "a"  "b"  
## t61  "b"  "a"  "b"  "a"  "b"  "a"  "a"  "b"  "a"  "b"  
## t62  "b"  "a"  "b"  "b"  "a"  "a"  "a"  "a"  "a"  "b"  
## t10  "b"  "a"  "b"  "b"  "b"  "b"  "a"  "b"  "a"  "b"  
## t11  "b"  "a"  "a"  "b"  "b"  "a"  "b"  "b"  "a"  "b"  
## t4   "b"  "b"  "b"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  
## t13  "b"  "a"  "a"  "a"  "b"  "a"  "a"  "a"  "b"  "b"  
## t31  "b"  "a"  "a"  "b"  "a"  "b"  "b"  "a"  "a"  "a"  
## t93  "b"  "a"  "a"  "b"  "b"  "b"  "b"  "a"  "b"  "b"  
## t94  "b"  "a"  "a"  "a"  "b"  "a"  "a"  "a"  "b"  "a"  
## t8   "b"  "a"  "a"  "b"  "b"  "b"  "a"  "a"  "a"  "b"  
## t5   "a"  "a"  "b"  "b"  "b"  "b"  "b"  "b"  "a"  "b"  
## t42  "b"  "b"  "b"  "b"  "b"  "a"  "a"  "a"  "a"  "a"  
## t76  "a"  "b"  "b"  "a"  "b"  "a"  "b"  "a"  "a"  "a"  
## t77  "a"  "b"  "b"  "b"  "a"  "b"  "b"  "b"  "a"  "b"  
## t63  "a"  "a"  "b"  "b"  "b"  "b"  "a"  "a"  "b"  "a"  
## t95  "b"  "b"  "b"  "b"  "a"  "b"  "a"  "a"  "a"  "a"  
## t96  "b"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  "b"  "a"  
## t32  "b"  "b"  "b"  "a"  "a"  "b"  "a"  "b"  "b"  "a"  
## t33  "b"  "b"  "a"  "b"  "b"  "b"  "b"  "b"  "b"  "a"  
## t81  "b"  "a"  "a"  "b"  "b"  "b"  "a"  "a"  "a"  "b"  
## t82  "b"  "a"  "a"  "b"  "b"  "b"  "a"  "a"  "a"  "b"  
## t24  "b"  "b"  "a"  "a"  "b"  "b"  "a"  "b"  "a"  "b"  
## t14  "a"  "a"  "a"  "a"  "a"  "a"  "a"  "b"  "a"  "b"  
## t25  "a"  "a"  "b"  "a"  "a"  "a"  "b"  "a"  "b"  "b"  
## t78  "a"  "a"  "b"  "b"  "a"  "a"  "b"  "b"  "b"  "b"  
## t79  "a"  "b"  "b"  "a"  "a"  "a"  "b"  "a"  "b"  "b"  
## t19  "b"  "b"  "a"  "b"  "b"  "a"  "a"  "b"  "a"  "a"  
## t50  "a"  "b"  "a"  "a"  "b"  "b"  "a"  "a"  "a"  "b"  
## t51  "b"  "a"  "b"  "a"  "a"  "a"  "a"  "a"  "a"  "a"  
## t26  "a"  "b"  "a"  "b"  "b"  "b"  "a"  "a"  "a"  "a"  
## t66  "a"  "a"  "b"  "a"  "b"  "a"  "a"  "b"  "a"  "a"  
## t67  "a"  "a"  "b"  "a"  "a"  "a"  "b"  "a"  "a"  "b"  
## t54  "a"  "b"  "b"  "a"  "b"  "a"  "b"  "a"  "a"  "a"  
## t55  "b"  "a"  "a"  "a"  "a"  "a"  "a"  "a"  "a"  "b"  
## t52  "a"  "a"  "a"  "a"  "b"  "b"  "b"  "b"  "b"  "b"  
## t53  "a"  "a"  "a"  "a"  "a"  "b"  "a"  "b"  "b"  "b"  
## t7   "b"  "a"  "b"  "a"  "b"  "a"  "a"  "a"  "b"  "b"  
## t18  "b"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  "b"  "b"  
## t20  "b"  "b"  "b"  "b"  "a"  "b"  "b"  "b"  "a"  "a"  
## t69  "b"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  "b"  "b"  
## t87  "b"  "b"  "b"  "b"  "b"  "b"  "a"  "b"  "b"  "b"  
## t88  "b"  "b"  "b"  "b"  "b"  "a"  "b"  "b"  "b"  "b"  
## t65  "b"  "b"  "b"  "b"  "a"  "b"  "b"  "b"  "a"  "a"  
## t74  "b"  "b"  "b"  "a"  "b"  "a"  "a"  "b"  "b"  "a"  
## t75  "a"  "a"  "a"  "a"  "a"  "b"  "a"  "b"  "a"  "b"  
## t39  "b"  "a"  "a"  "a"  "b"  "b"  "b"  "a"  "a"  "b"  
## t22  "a"  "a"  "a"  "a"  "a"  "b"  "b"  "a"  "a"  "b"  
## t34  "b"  "b"  "b"  "a"  "b"  "b"  "a"  "b"  "a"  "a"  
## t35  "a"  "a"  "a"  "b"  "b"  "a"  "b"  "a"  "a"  "b"  
## t70  "b"  "a"  "a"  "a"  "b"  "b"  "b"  "b"  "b"  "a"  
## t71  "a"  "a"  "b"  "a"  "b"  "a"  "a"  "a"  "a"  "b"  
## t99  "a"  "b"  "b"  "a"  "b"  "a"  "a"  "b"  "b"  "b"  
## t100 "b"  "a"  "a"  "a"  "b"  "a"  "a"  "b"  "a"  "b"  
## t80  "b"  "a"  "a"  "a"  "b"  "a"  "a"  "b"  "a"  "b"  
## t6   "a"  "a"  "a"  "a"  "a"  "b"  "a"  "b"  "a"  "a"  
## t12  "b"  "a"  "a"  "a"  "a"  "b"  "a"  "a"  "b"  "b"  
## t30  "b"  "a"  "a"  "b"  "a"  "a"  "a"  "b"  "b"  "b"  
## t56  "b"  "b"  "a"  "a"  "b"  "b"  "a"  "a"  "b"  "b"  
## t57  "b"  "a"  "a"  "b"  "a"  "b"  "a"  "b"  "b"  "a"  
## t40  "a"  "b"  "b"  "a"  "a"  "b"  "b"  "b"  "a"  "b"  
## t41  "b"  "a"  "b"  "a"  "b"  "b"  "a"  "b"  "b"  "a"  
## t1   "a"  "b"  "b"  "a"  "b"  "a"  "a"  "b"  "b"  "b"  
## t9   "a"  "b"  "b"  "a"  "a"  "b"  "a"  "b"  "b"  "b"  
## t48  "b"  "a"  "a"  "a"  "b"  "b"  "a"  "b"  "b"  "b"  
## t49  "b"  "a"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  "b"  
## t72  "b"  "a"  "a"  "a"  "a"  "b"  "a"  "b"  "b"  "b"  
## t73  "a"  "a"  "a"  "b"  "b"  "b"  "a"  "b"  "a"  "b"

Note that it is arbitrary, and unnecessary, that I have simulated using the same transition matrix Q for each data column here.

Finally, our tree-transformation fitting function. Here, I've used ace from ape internally to compute the likelihood of each character for each tree transformation. One could instead use fitDiscrete (which is slower, but perhaps more robust) or function in the phangorn package.

## here is our model fitting function
fitTransform<-function(tree,X,model,interval){
    lk<-function(par,tree,X,model){
        sum(apply(X,2,function(x,tree,model,par) 
            logLik(ace(x,rescale(tree,model,par),type="discrete")),
            tree=tree,model=model,par=par))
    }
    fit<-optimize(lk,interval,tree=tree,X=X,model=model,maximum=TRUE)
    return(list(par=fit$maximum,model=model,logLik=fit$objective))
}
obj<-fitTransform(tree,X,model="lambda",interval=c(0,1))
## here is our fitted model
obj

## $par
## [1] 0.472
## 
## $model
## [1] "lambda"
## 
## $logLik
## [1] -666.7

We could use the same function with different values of model and (necessarily, if so) interval. I use model="lambda" merely to illustrate the case.

That's it.