T-shirt design for Ilhabela workshop - in R

I'm in the midst of working on a design for the t-shirts for this year's macroevolution workshop in Ilhabela, Brazil. I've created a tentative 'front' design entirely in R using the function of phytools, as follows:

library(phytools)
## don't run
# tree<-pbtree(n=100)
# x<-fastBM(tree)

obj<-contMap(tree,x,plot=FALSE)
obj<-setMap(obj,colors=c("blue","purple","red"))
layout(mat=matrix(c(1,2),2,1),heights=c(0.8,0.2))
par(bg="black")
par(fg="white")
plotSimmap(paintSubTree(tree,Ntip(tree)+1,"1"),type="fan",
    ftype="off",colors=setNames("white","1"),lwd=6,part=0.5)

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

plotSimmap(obj$tree,type="fan",ftype="off",colors=obj$cols,lwd=4,
    add=TRUE,part=0.5)

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

plot.new()
text(0.5,0.5,"Latin American Macroevolution Workshop\nIlhabela Brazil 2015",
    col="white",cex=2.1,font=2)

The knitr output is kind of aliased. We can create a higher quality PDF easily:

pdf(file="t-shirt.pdf",width=8,height=4.75)
layout(mat=matrix(c(1,2),2,1),heights=c(0.8,0.2))
par(bg="black")
par(fg="white")
plotSimmap(paintSubTree(tree,Ntip(tree)+1,"1"),type="fan",
    ftype="off",colors=setNames("white","1"),lwd=6,part=0.5)

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

plotSimmap(obj$tree,type="fan",ftype="off",colors=obj$cols,lwd=4,
    add=TRUE,part=0.5)

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

plot.new()
text(0.5,0.5,"Latin American Macroevolution Workshop\nIlhabela Brazil 2015",
    col="white",cex=2.1,font=2)
dev.off()

## windows 
##       2

Here is a link to the file. (It looks much better!)

That's it for now.

Testing a null hypothesis of K=1.0 with measurement error in the estimation of trait means for species

In the past I've described a “parametric bootstrapping” approach for testing the null hypothesis that Blomberg's K is significantly different from 1.0 (for instance here). 1.0 is the expected value of Blomberg's K under a Brownian motion model of evolutionary change for our character. Here, I'll describe the relatively simple modification of that procedure that can be used to conduct the same hypothesis test under conditions in which we also have sampling error in the estimation of species means for our character.

Let's start by simulating some trait data with sampling error. We can do this using the phytools functions pbtree (for phylogeny simulation), fastBM (to simulate trait data), and sampleFrom, which is a simple function, normally used internally, which allows us to sample trait data under some model (for instance, with error).

library(phytools)
## simulate tree
tree<-pbtree(n=100,scale=1)
## simulate trait data
x<-fastBM(tree,sig2=1)
## add error
ve<-setNames(rchisq(n=Ntip(tree),df=1)/runif(n=Ntip(tree),min=1,max=5),
    tree$tip.label)
xe<-sampleFrom(x,ve,n=rep(1,Ntip(tree)))
plot(x,xe,xlab="true species means",
    ylab="simulated species means with sampling error")

Now let's fit our model & generate data under the null hypothesis of K = 1.0. This is where we would normally begin with a genuine empirical dataset.

fit<-phylosig(tree,xe,se=sqrt(ve))
fit

## $K
## [1] 0.9085662
## 
## $sig2
## [1] 1.125822
## 
## $logL
## [1] -112.5553

nullX<-fastBM(tree,nsim=200,sig2=fit$sig2)
nullXe<-apply(nullX,2,sampleFrom,xvar=ve,n=rep(1,Ntip(tree)))
obj<-apply(nullXe,2,phylosig,tree=tree,se=sqrt(ve))
nullK<-sapply(obj,function(x) x[[1]])
hist(nullK,breaks=20,xlab="null distribution for K",
    main="Null distribution of K")
lines(c(fit$K,fit$K),c(0,par()$usr[4]),lty="dashed",col="red")
text(x=fit$K,y=0.985*par()$usr[4],"observed value of K",
    pos=4,offset=0.2)

Obviously, here we are well within the null distribution for K, but we can also attach a P-value to this observation. I usually use the logarithm of K, because this equally penalizes K = 0.5 and K = 2.0.

P<-mean(abs(log(nullK))>=abs(log(fit$K))) ## two-tailed test
P

## [1] 0.675

That's it.

Simple method to rescale tree to a particular mean tip height

Here's a simple trick to rescale a tree to have a particular mean height (if not ultrametric):

library(phytools)
## in this case we will use a random tree
tree<-rtree(n=26,tip.label=LETTERS)
plotTree(tree,mar=c(5.1,0.1,0.1,0.1))
axis(1)
## compute all tip heights
h<-sapply(1:Ntip(tree),nodeheight,tree=tree)
h

##  [1] 2.998386 3.391986 3.502868 3.567715 3.110754 3.365106 2.714611
##  [8] 1.358671 2.152025 2.313791 2.715758 3.055251 3.199052 1.739875
## [15] 3.869145 4.121823 3.534816 3.962804 1.561411 2.536063 3.378451
## [22] 2.719116 2.628788 2.661650 2.900916 2.311506

mean(h)

## [1] 2.898936

lines(c(mean(h),mean(h)),par()$usr[3:4],lty="dashed",col="red")

## decide on new desired mean height
new.h<-100
## rescale tree
tree$edge.length<-tree$edge.length/mean(h)*new.h
h<-sapply(1:Ntip(tree),nodeheight,tree=tree)
h

##  [1] 103.43058 117.00797 120.83289 123.06981 107.30675 116.08072  93.64164
##  [8]  46.86791  74.23498  79.81517  93.68118 105.39215 110.35264  60.01771
## [15] 133.46776 142.18398 121.93495 136.69855  53.86151  87.48253 116.54105
## [22]  93.79704  90.68114  91.81473 100.06830  79.73635

mean(h)

## [1] 100

plotTree(tree,mar=c(5.1,0.1,0.1,0.1))
axis(1)
lines(c(mean(h),mean(h)),par()$usr[3:4],lty="dashed",col="red")

That's it.