Analysis for Revell et al. (2018; Methods in Ecol. & Evol.) Analysis for:
Liam J. Revell, Laura E. González-Valenzuela, Alejandro Alfonso, Luisa A. Castellanos-García, Carlos E. Guarnizo, and Andrew J. Crawford. 2018. Comparing evolutionary rates between trees, clades, & traits. Methods in Ecology & Evolution.

Preamble:

The following code reproduces exactly the analysis of Revell et al. (2018; MEE) using R version:

R.version
##                _                           
## platform       x86_64-w64-mingw32          
## arch           x86_64                      
## os             mingw32                     
## system         x86_64, mingw32             
## status                                     
## major          3                           
## minor          4.3                         
## year           2017                        
## month          11                          
## day            30                          
## svn rev        73796                       
## language       R                           
## version.string R version 3.4.3 (2017-11-30)
## nickname       Kite-Eating Tree

and phytools version:

packageVersion("phytools")
## [1] '0.6.53'

although it will probably work with earlier versions of phytools & R.

Figure 1:

Figure 1 & Table 1 give the results from a type I error analysis of the method, using both a χ2 distribution for the test statistic & a distribution obtained via simulation under the null hypothesis. The simulations & results follow.

Simulations:

library(phytools)
## Loading required package: ape
## Loading required package: maps
set.seed(1)

nrep<-1000
t10<-replicate(nrep,pbtree(n=10,nsim=2),simplify=FALSE)
x10<-lapply(t10,function(x) lapply(x,fastBM))
t50<-replicate(nrep,pbtree(n=50,nsim=2),simplify=FALSE)
x50<-lapply(t50,function(x) lapply(x,fastBM))
foo<-function(n1,n2) c(pbtree(n=n1),pbtree(n=n2))
t10.t50<-replicate(nrep,foo(10,50),simplify=FALSE)
x10.x50<-lapply(t10.t50,function(x) lapply(x,fastBM))
t30<-replicate(nrep,pbtree(n=30,nsim=3),simplify=FALSE)
x30<-lapply(t30,function(x) lapply(x,fastBM))

typeIerror<-function(i,t,x,max){
    obj<-ratebytree(t,x)
    if(i==1) cat("|.") else cat(".")
    if(i==max) cat("|\n")
    if(i%%50==0) cat("\n ")
    flush.console()
    obj$P.chisq
}

Results:

par(mfrow=c(2,2))

P10<-mapply(typeIerror,1:nrep,t10,x10,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
h10<-hist(P10,breaks=seq(0,1,by=0.05),col="grey",xlab="P-value",
    main="",freq=FALSE)
mtext(text=expression('a) N'[1]*'=N'[2]*'=10'),adj=0,line=1,
    cex=1)
abline(h=1,lwd=2,col="red",lty="dashed")

P50<-mapply(typeIerror,1:nrep,t50,x50,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
h50<-hist(P50,breaks=seq(0,1,by=0.05),col="grey",main="",
    xlab="P-value",freq=FALSE)
mtext(text=expression('b) N'[1]*'=N'[2]*'=50'),adj=0,line=1,
    cex=1)
abline(h=1,lwd=2,col="red",lty="dashed")

P10.50<-mapply(typeIerror,1:nrep,t10.t50,x10.x50,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
h10.50<-hist(P10.50,breaks=seq(0,1,by=0.05),col="grey",main="",freq=FALSE,
    xlab="P-value")
mtext(text=expression('c) N'[1]*'=10, N'[2]*'=50'),adj=0,line=1,
    cex=1)
abline(h=1,lwd=2,col="red",lty="dashed")

P30<-mapply(typeIerror,1:nrep,t30,x30,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
h30<-hist(P30,breaks=seq(0,1,by=0.05),col="grey",freq=FALSE,main="",
    xlab="P-value")
mtext(text=expression('d) N'[1]*'=N'[2]*'=N'[3]*'=30'),
    adj=0,line=1,cex=1)
abline(h=1,lwd=2,col="red",lty="dashed")

plot of chunk unnamed-chunk-4

Figure 1. Distribution of P-values obtained from hypothesis tests when data were simulated under the null hypothesis of no difference in rate between trees. The expected distribution is uniform on the interval of [0, 1]. Panels a) through d) show the number of species in the sets of two or three trees simulated for each condition.

Table 1:

Table 1 summarizes the type I error from both the parametric test and the test conducted via simulation under the null hypothesis.

typeI<-matrix(sapply(list(P10,P50,P10.50,P30),function(x,alpha)
    mean(x<=alpha),alpha=0.05),4,1,dimnames=list(
    c("N1=N2=10","N1=N2=50","N1=10,N2=50","N1=N2=N3=30"),
    "type I error"))
typeI<-cbind(typeI,sapply(typeI[,1]*nrep,pbinom,size=nrep,prob=0.05,
    lower.tail=FALSE))
colnames(typeI)[2]<-"P(binomial)"

## now repeat type I error test with simulation
typeIerror.sim<-function(i,t,x,max){
    obj<-ratebytree(t,x,test="simulation",quiet=TRUE)
    if(i==1) cat("|.") else cat(".")
    if(i==max) cat("|\n")
    if(i%%50==0) cat("\n ")
    flush.console()
    obj$P.sim
}

P10.sim<-mapply(typeIerror.sim,1:nrep,t10,x10,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
P50.sim<-mapply(typeIerror.sim,1:nrep,t50,x50,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  .......................
## Warning in sqrt(diag(Cov.multi)[1:N + m]): NaNs produced
## ...........................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
P10.50.sim<-mapply(typeIerror.sim,1:nrep,t10.t50,x10.x50,
    MoreArgs=list(max=nrep))
## |..................................................
##  ....................................
## Warning in sqrt(diag(Cov.multi)[1:N + m]): NaNs produced
## Warning in sqrt(diag(Cov.onerate)[1:N + 1]): NaNs produced
## ..............
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
P30.sim<-mapply(typeIerror.sim,1:nrep,t30,x30,MoreArgs=list(max=nrep))
## |..................................................
##  ..................................................
##  ..................................................
##  ....
## Warning in sqrt(diag(Cov.multi)[1:N + m]): NaNs produced

## Warning in sqrt(diag(Cov.multi)[1:N + m]): NaNs produced
## ..............................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................
##  ..................................................|
## 
## 
typeI.sim<-matrix(sapply(list(P10.sim,P50.sim,P10.50.sim,P30.sim),
    function(x,alpha) mean(x<=alpha),alpha=0.05),4,1,
    dimnames=list(
    c("N1=N2=10","N1=N2=50","N1=10,N2=50","N1=N2=N3=30"),
    "type I error"))
typeI.sim<-cbind(typeI.sim,sapply(typeI.sim[,1]*nrep,pbinom,size=nrep,
    prob=0.05,lower.tail=FALSE))
colnames(typeI.sim)[2]<-"P(binomial)"

typeI<-cbind(typeI,typeI.sim)
colnames(typeI)<-c("type I (X^2)","P(binomial)",
    "type I (simulation)","P(binomial)")
print(typeI,digits=4)
##             type I (X^2) P(binomial) type I (simulation) P(binomial)
## N1=N2=10           0.065     0.01493               0.045      0.7390
## N1=N2=50           0.052     0.35140               0.052      0.3514
## N1=10,N2=50        0.056     0.17206               0.043      0.8267
## N1=N2=N3=30        0.063     0.02843               0.047      0.6344

Table 1. Type I error rates of a likelihood method of comparing evolutionary rates between trees when the null distribution for the likelihood-ratio test statistic was either: a χ2 distribution with degrees of freedom equal to the number of trees, k, minus 1; or a distribution obtained via simulation under the null hypothesis of no difference in rate. P-values were obtained by comparison of the observed type I error rate to a binomial distribution with probability set to the level of α employed in the statistical test (0.05).

Figure 2:

Figure 2 gives a projection of the tree into morphospace of three simulated phylogenies: two sharing a common rate of evolution, and the third with a different rate.

set.seed(10)

t1<-pbtree(n=30)
t2<-pbtree(n=30)
t3<-pbtree(n=30)
x<-fastBM(t1,sig2=1)
y<-fastBM(t2,sig2=2)
z<-fastBM(t3,sig2=1)
trees<-c(t1,t2,t3)
x<-list(x,y,z)

ylim<-range(x)
par(mfrow=c(1,3))
phenogram(trees[[1]],x[[1]],ylim=ylim,spread.cost=c(1,0))
mtext(text=expression('a) '*sigma[1]^2*' = 1.0'),adj=0,line=1,cex=1)
phenogram(trees[[2]],x[[2]],ylim=ylim,spread.cost=c(1,0))
mtext(text=expression('b) '*sigma[2]^2*' = 2.0'),adj=0,line=1,cex=1)
phenogram(trees[[3]],x[[3]],ylim=ylim,spread.cost=c(1,0))
mtext(text=expression('c) '*sigma[3]^2*' = 1.0'),adj=0,line=1,cex=1)

plot of chunk unnamed-chunk-6

Figure 2. Example simulated trees and data for N1=N2=N3=30 species in which Brownian rate parameters, σ1232=1.0 and σ22=2.0. Each panel shows a traitgram, i.e., a projection of the phylogeny into phenotypic space. Note that the trees were simulated under a scenario of pure-birth with a constant rate of speciation and a taxon-number stopping condition, thus resulting in different total depths.

Figures 3:

Figure 3 gives the power of the method for various genuine differences in the simulated rate of evolution between trees. Simulations & results follow.

Simulations:

set.seed(1)

nrep<-100
t10<-replicate(nrep,pbtree(n=10,nsim=2),simplify=FALSE)
t50<-replicate(nrep,pbtree(n=50,nsim=2),simplify=FALSE)
foo<-function(n1,n2) c(pbtree(n=n1),pbtree(n=n2))
t10.t50<-replicate(nrep,foo(10,50),simplify=FALSE)
t30<-replicate(nrep,pbtree(n=30,nsim=3),simplify=FALSE)

power<-function(t,sig2){
    x<-fastBM(t[[1]])
    y<-fastBM(t[[2]],sig2=sig2)
    z<-if(length(t)==3) fastBM(t[[3]]) else NULL
    x<-if(is.null(z)) list(x,y) else list(x,y,z)
    obj<-ratebytree(t,x)
    if(is.null(z)) 
        setNames(c(obj$multi.rate.model$sig2,
            obj$multi.rate.model$logL,
            obj$P.chisq),
            c("sig2[1]","sig2[2]","log(L)","P"))
    else 
        setNames(c(obj$multi.rate.model$sig2,
            obj$multi.rate.model$logL,
            obj$P.chisq),
            c("sig2[1]","sig2[2]","sig2[3]","log(L)",
            "P"))
}
foo<-function(sig2,t) t(sapply(t,power,sig2=sig2))

sig2<-seq(1,4,by=0.2)
fitted.t10<-setNames(lapply(sig2,foo,t=t10),sig2)
fitted.t50<-setNames(lapply(sig2,foo,t=t50),sig2)
fitted.t10.t50<-setNames(lapply(sig2,foo,t=t10.t50),sig2)
fitted.t30<-setNames(lapply(sig2,foo,t=t30),sig2)
par(mfrow=c(2,2))

plot(sig2,sapply(fitted.t10,function(x) mean(x[,"P"]<=0.05)),
    type="b",xlab=expression(sigma[2]^2/sigma[1]^2),ylab="Power",
    ylim=c(0,1),lwd=2)
mtext(text=expression('a) N'[1]*'=N'[2]*'=10'),adj=0,line=1,
    cex=1)
abline(h=0.05,lty="dashed")
abline(h=1,lty="dashed")
text(x=3.9,y=0.05,"0.05",pos=3)
text(x=1.1,y=1,"1.00",pos=1)

plot(sig2,sapply(fitted.t50,function(x) mean(x[,"P"]<=0.05)),
    type="b",xlab=expression(sigma[2]^2/sigma[1]^2),ylab="Power",
    ylim=c(0,1),lwd=2)
mtext(text=expression('b) N'[1]*'=N'[2]*'=50'),adj=0,line=1,
    cex=1)
abline(h=0.05,lty="dashed")
abline(h=1,lty="dashed")
text(x=3.9,y=0.05,"0.05",pos=3)
text(x=1.1,y=1,"1.00",pos=1)

plot(sig2,sapply(fitted.t10.t50,function(x) mean(x[,"P"]<=0.05)),
    type="b",xlab=expression(sigma[2]^2/sigma[1]^2),ylab="Power",
    ylim=c(0,1),lwd=2)
mtext(text=expression('c) N'[1]*'=10, N'[2]*'=50'),adj=0,line=1,
    cex=1)
abline(h=0.05,lty="dashed")
abline(h=1,lty="dashed")
text(x=3.9,y=0.05,"0.05",pos=3)
text(x=1.1,y=1,"1.00",pos=1)

plot(sig2,sapply(fitted.t30,function(x) mean(x[,"P"]<=0.05)),
    type="b",xlab=expression(sigma[2]^2/sigma[1]^2),ylab="Power",
    ylim=c(0,1),lwd=2)
mtext(text=expression('d) N'[1]*'=N'[2]*'=N'[3]*'=30'),
    adj=0,line=1,cex=1)
abline(h=0.05,lty="dashed")
abline(h=1,lty="dashed")
text(x=3.9,y=0.05,"0.05",pos=3)
text(x=1.1,y=1,"1.00",pos=1)

plot of chunk unnamed-chunk-8

Figure 3. Power to reject the null hypothesis of equal rates under different simulation conditions of the Brownian rate parameters, σ12 and σ22 (and σ32, for panel d). Invariably, the generating value of σ12 (and σ32, if applicable) was fixed at 1.0, while σ22 was varied between 1.0 and 4.0 in intervals of 0.2. Panels a) through d) are as in Figure 1.

Figure 4:

Figure 4 re-uses the simulations from Figure 3, but to conduct an analysis of parameter estimation when trees genuinely differ in their rate of evolution.

library(plotrix)

par(mfrow=c(2,2))
par(mar=c(5.1,5.1,4.1,1.1))

sd1<-sapply(fitted.t10,function(x) sqrt(var(x[,"sig2[1]"])))
sd2<-sapply(fitted.t10,function(x) sqrt(var(x[,"sig2[2]"])))
plotCI(sig2,sapply(fitted.t10,function(x) mean(x[,"sig2[2]"])),
    uiw=sd2,liw=sd2,ylab=expression('estimated '*sigma[1]^2*' & '*
    sigma[2]^2),xlab=expression('simulated '*sigma[2]^2/sigma[1]^2),
    lwd=1,gap=0.04,sfrac=0)
points(sig2,sapply(fitted.t10,function(x) mean(x[,"sig2[2]"])),
    pch=21,bg="grey")
plotCI(sig2,sapply(fitted.t10,function(x) mean(x[,"sig2[1]"])),
    uiw=sd1,liw=sd1,lwd=1,gap=0.04,add=TRUE,sfrac=0)
points(sig2,sapply(fitted.t10,function(x) mean(x[,"sig2[1]"])),
    pch=21,bg="grey")
mtext(text=expression('a) N'[1]*'=N'[2]*'=10'),adj=0,line=1,
    cex=1)
lines(c(0,5),c(0,5),col="red",lwd=2,lty="dashed")
lines(c(0,5),c(1,1),col="blue",lwd=2,lty="dashed")

sd1<-sapply(fitted.t50,function(x) sqrt(var(x[,"sig2[1]"])))
sd2<-sapply(fitted.t50,function(x) sqrt(var(x[,"sig2[2]"])))
plotCI(sig2,sapply(fitted.t50,function(x) mean(x[,"sig2[2]"])),
    uiw=sd2,liw=sd2,ylab=expression('estimated '*sigma[1]^2*' & '*
    sigma[2]^2),xlab=expression('simulated '*sigma[2]^2/sigma[1]^2),
    lwd=1,gap=0.04,sfrac=0)
points(sig2,sapply(fitted.t50,function(x) mean(x[,"sig2[2]"])),
    pch=21,bg="grey")
plotCI(sig2,sapply(fitted.t50,function(x) mean(x[,"sig2[1]"])),
    uiw=sd1,liw=sd1,lwd=1,gap=0.04,add=TRUE,sfrac=0)
points(sig2,sapply(fitted.t50,function(x) mean(x[,"sig2[1]"])),
    pch=21,bg="grey")
mtext(text=expression('b) N'[1]*'=N'[2]*'=50'),adj=0,line=1,
    cex=1)
lines(c(0,5),c(0,5),col="red",lwd=2,lty="dashed")
lines(c(0,5),c(1,1),col="blue",lwd=2,lty="dashed")

sd1<-sapply(fitted.t10.t50,function(x) sqrt(var(x[,"sig2[1]"])))
sd2<-sapply(fitted.t10.t50,function(x) sqrt(var(x[,"sig2[2]"])))
plotCI(sig2,sapply(fitted.t10.t50,function(x) mean(x[,"sig2[2]"])),
    uiw=sd2,liw=sd2,ylab=expression('estimated '*sigma[1]^2*' & '*
    sigma[2]^2),xlab=expression('simulated '*sigma[2]^2/sigma[1]^2),
    lwd=1,gap=0.04,sfrac=0)
points(sig2,sapply(fitted.t10.t50,function(x) mean(x[,"sig2[2]"])),
    pch=21,bg="grey")
plotCI(sig2,sapply(fitted.t10.t50,function(x) mean(x[,"sig2[1]"])),
    uiw=sd1,liw=sd1,lwd=1,gap=0.04,add=TRUE,sfrac=0)
points(sig2,sapply(fitted.t10.t50,function(x) mean(x[,"sig2[1]"])),
    pch=21,bg="grey")
mtext(text=expression('c) N'[1]*'=10, N'[2]*'=50'),adj=0,line=1,
    cex=1)
lines(c(0,5),c(0,5),col="red",lwd=2,lty="dashed")
lines(c(0,5),c(1,1),col="blue",lwd=2,lty="dashed")

sd1<-sapply(fitted.t30,function(x) sqrt(var(x[,"sig2[1]"])))
sd2<-sapply(fitted.t30,function(x) sqrt(var(x[,"sig2[2]"])))
sd3<-sapply(fitted.t30,function(x) sqrt(var(x[,"sig2[3]"])))
plotCI(sig2,sapply(fitted.t30,function(x) mean(x[,"sig2[2]"])),
    uiw=sd2,liw=sd2,ylab=expression('estimated '*sigma[1]^2*', '*
    sigma[2]^2*', & '*sigma[3]^2),
    xlab=expression('simulated '*sigma[2]^2/sigma[1]^2),
    lwd=1,gap=0.04,sfrac=0)
points(sig2,sapply(fitted.t30,function(x) mean(x[,"sig2[2]"])),
    pch=21,bg="grey")
plotCI(sig2,sapply(fitted.t30,function(x) mean(x[,"sig2[1]"])),
    uiw=sd1,liw=sd1,lwd=1,gap=0.04,add=TRUE,sfrac=0)
points(sig2,sapply(fitted.t30,function(x) mean(x[,"sig2[1]"])),
    pch=21,bg="grey")
plotCI(sig2,sapply(fitted.t30,function(x) mean(x[,"sig2[3]"])),
    uiw=sd1,liw=sd1,lwd=1,gap=0.04,add=TRUE,sfrac=0)
points(sig2,sapply(fitted.t30,function(x) mean(x[,"sig2[3]"])),
    pch=21,bg="grey")
mtext(text=expression('d) N'[1]*'=N'[2]*'=N'[3]*'=30'),
    adj=0,line=1,cex=1)
lines(c(0,5),c(0,5),col="red",lwd=2,lty="dashed")
lines(c(0,5),c(1,1),col="blue",lwd=2,lty="dashed")

plot of chunk unnamed-chunk-9

Figure 4. Parameter estimation for simulation conditions involving a genuine difference in rate between trees. As in Figure 2, the generating value of the Brownian rate parameter σ12 (and σ32, if applicable) was fixed at 1.0 while σ22 varied between 1.0 and 4.0. For all panels the blue line represents the generating values of σ12 (and σ32, if applicable) and the red line represents the generating values of σ22. Panels a) through d) are as in Figures 1 and 2.

Figure 5:

Figure 5 shows a projection of the phylogeny into a space defined by body size on the vertical axis and time since the root on the horizontal axis for two different lizard clades. We will then analyze these data, and data for body shape, for a difference in rate between clades. In the following code, we read data from file, visualize the data for body size evolution, extra PC 1 from a phylogenetic PCA, and analyze both traits for a difference in rate between trees. The data are from Harmon et al. (2010).

library(geiger)
## 
## Attaching package: 'geiger'
## The following object is masked from 'package:plotrix':
## 
##     rescale
ages<-setNames(c(75,74),c("phryn","liol"))
ages
## phryn  liol 
##    75    74
## a simple function to rescale
rescaleTree<-function(t,d){
    h<-max(nodeHeights(t))
    t$edge.length<-t$edge.length/h*d
    t
}

Analysis:

phryn.tree<-read.tree("phrynosomatpl.phy")
phryn.tree<-rescaleTree(phryn.tree,ages["phryn"]/100)
phryn.data<-read.csv("iguanids.csv",row.names=1,header=TRUE,
    na.string=".")
chk<-name.check(phryn.tree,phryn.data)
phryn.data<-phryn.data[phryn.tree$tip.label,]
phryn.size<-log(as.matrix(phryn.data)[,1]) ## size on a log scale.
liol.tree<-read.tree("liolaeminipl.phy")
liol.tree<-rescaleTree(liol.tree,ages["liol"]/100)
liol.data<-read.csv("iguanids.csv",row.names=1,header=TRUE,
    na.string=".")
chk<-name.check(liol.tree,liol.data)
liol.data<-liol.data[liol.tree$tip.label,]
liol.size<-log(as.matrix(liol.data)[,1])
a.phryn<-exp(fastAnc(phryn.tree,phryn.size))
a.liol<-exp(fastAnc(liol.tree,liol.size))
par(mfrow=c(1,2),mar=c(5.1,4.1,3.1,1.1))
phenogram(phryn.tree,c(exp(phryn.size),a.phryn),
    spread.cost=c(1,0),ftype="i",fsize=0.6,
    ylab="SVL(mm)",xlab="time(ma x 100)",col="lightblue",
    ylim=range(exp(c(phryn.size,liol.size))))
## Optimizing the positions of the tip labels...
mtext(text="a) Phrynosomatinae body size",cex=1.2,at=0.3)
labels<-setNames(seq(min(exp(c(phryn.size,liol.size))),
    max(exp(c(phryn.size,liol.size))),
    by=diff(range(exp(c(phryn.size,liol.size))))/(Ntip(liol.tree)-1)),
    names(sort(exp(liol.size))))
phenogram(liol.tree,c(exp(liol.size),a.liol),
    spread.cost=c(1,0),ftype="i",fsize=0.6,
    ylab="SVL(mm)",xlab="time(ma x 100)",col="lightgreen",
    ylim=range(exp(c(phryn.size,liol.size))),
    label.pos=labels)
## Optimizing the positions of the tip labels...
mtext(text="b) Liolaemini body size",cex=1.2,at=0.2)

plot of chunk unnamed-chunk-11

Figure 5. Projection of the tree into phenotype space (traitgram) for body size evolution in two lizard clades: the North American phrynosomatines; and the South American tribe, Liolaemini. Rate estimates and ancestral states for the projection were obtained on a natural logarithm scale and then back translated onto a linear scale for the purposes of visualization only.

Table 2

Table 2, here shown in two parts, summarizes the results of a comparison in rate between trees for body size & body shape, as defined above.

fit.size<-ratebytree(c(phryn.tree,liol.tree),list(phryn.size,
    liol.size))
fit.size
## ML common-rate model:
##  s^2  a[1]   a[2]    k   logL
## value    0.2613  4.1772  4.2621  3   -4.8496 
## SE       0.0316  0.1478  0.2394  
## 
## ML multi-rate model:
##   s^2[1] s^2[2]   a[1]   a[2]    k   logL
## value    0.1911  0.3346  4.1772  4.2621  4   -2.1877 
## SE       0.0323  0.0578  0.1264  0.2709  
## 
## Likelihood ratio: 5.3238 
## P-value (based on X^2): 0.021 
## 
## R thinks it has found the ML solution.

Table 2.1. Body size (SVL) for two lizard clades: (1) the North American iguanian subfamily Phrynosomatinae; and (2) the South American lizard tribe Liolaemini.

taxa.na<-names(which(apply(phryn.data,1,function(x) any(is.na(x)))))
phryn.tree<-drop.tip(phryn.tree,taxa.na)
phryn.data<-phryn.data[phryn.tree$tip.label,]
pca.phryn<-phyl.pca(phryn.tree,log(phryn.data))
taxa.na<-names(which(apply(liol.data,1,function(x) any(is.na(x)))))
liol.tree<-drop.tip(liol.tree,taxa.na)
liol.data<-liol.data[liol.tree$tip.label,]
pca.liol<-phyl.pca(liol.tree,log(liol.data))
Evec<-(pca.liol$Evec+pca.phryn$Evec)/2
phryn.shape<-((as.matrix(log(phryn.data))-matrix(1,nrow(phryn.data),1)%*%
    t(phyl.vcv(as.matrix(log(phryn.data)),vcv(phryn.tree),1)$alpha))%*%
    Evec)[,2]
liol.shape<-((as.matrix(log(liol.data))-matrix(1,nrow(liol.data),1)%*%
    t(phyl.vcv(as.matrix(log(liol.data)),vcv(liol.tree),1)$alpha))%*%
    Evec)[,2]
fit.shape<-ratebytree(c(phryn.tree,liol.tree),list(phryn.shape,
    liol.shape))
fit.shape
## ML common-rate model:
##  s^2  a[1]   a[2]    k   logL
## value    0.1148  0   0   3   47.4324 
## SE       0.0144  0.0993  0.1588  
## 
## ML multi-rate model:
##   s^2[1] s^2[2]   a[1]   a[2]    k   logL
## value    0.1181  0.1112  0   0   4   47.4611 
## SE       0.0204  0.0201  0.1007  0.1563  
## 
## Likelihood ratio: 0.0574 
## P-value (based on X^2): 0.8106 
## 
## R thinks it has found the ML solution.

Table 2.2. Body shape (common phylogenetic PC2) for two lizard clades: (1) the North American iguanian subfamily Phrynosomatinae; and (2) the South American lizard tribe Liolaemini.