Computing phylogenetic signal from a set of quantitative traits in an efficient way

Today I received the following inquiry (simplified a bit below):

“I'm working with a large… dataset. I'm trying to find a way to apply the phylosig function across each column in my dataset so that I don't have to calculate K or λ one-by-one. Ideally I'd like to save the output value with the column header (the trait name). It seems like it should be pretty straightforward, but I have struggled to make this work, I think because the rownames must be assigned to each extracted data vector in order to match to the tree. I and some other R-users have several approaches (e.g., loops to extract individual vectors, convert to matrix and add rownames, various apply functions, etc). We haven't been successful.”

To show how to do this I'm going to use some real data from a study I published with Luke Mahler, Graham Reynolds, and Graham Slater a few years ago.

library(phytools)
Anolis.tree<-read.tree(file=
    "https://datadryad.org/bitstream/handle/10255/dryad.80249/Revell-etal.tree.tre")
Anolis.data<-read.csv(file=
    "https://datadryad.org/bitstream/handle/10255/dryad.80250/Revell-etal.data.csv",
    row.names=1)

We technically don't need to do this (as phylosig should take care of this for us), but just to avoid any error messages I'm going to 'clean up' my data frame so that the taxa in the data & the tree match exactly:

library(geiger)
chk<-name.check(Anolis.tree,Anolis.data)
chk

## $tree_not_data
## character(0)
## 
## $data_not_tree
## [1] "roosevelti"

Anolis.data<-Anolis.data[-which(rownames(Anolis.data)%in%
    chk$data_not_tree),]
name.check(Anolis.tree,Anolis.data)

## [1] "OK"

Just for fun, let's plot one of our characters: overall body size (SVL in mm).

plotTree.barplot(Anolis.tree,exp(Anolis.data[,1,drop=FALSE]),
    args.plotTree=list(ftype="off"),
    args.barplot=list(xlab="SVL (mm)",space=0.5))

Now all the characters together:

phylo.heatmap(Anolis.tree,Anolis.data,standardize=TRUE,
    fsize=c(0.5,0.7,1))

Ooh. That's cool. Note that the pattern looks pretty similar across all the traits because they are all highly correlated with body size!

Now let's compute our phylogenetic signals first using Blomberg et al.'s (2003) K:

K<-apply(Anolis.data,2,phylosig,tree=Anolis.tree)
K

##                AVG.SVL                 AVG.hl                 AVG.hw 
##               1.553678               1.609666               1.593687 
##                 AVG.hh                AVG.ljl           AVG.outlever 
##               1.633474               1.536071               1.565743 
## AVG.jugal.to.symphysis              AVG.femur              AVG.tibia 
##               1.570188               1.508540               1.520925 
##                AVG.met            AVG.ltoe.IV   AVG.toe.IV.lam.width 
##               1.585221               1.525271               1.450562 
##            AVG.humerus             AVG.radius           AVG.lfing.IV 
##               1.651481               1.647935               1.615006 
##  AVG.fing.IV.lam.width            AVG.pelv.ht            AVG.pelv.wd 
##               1.456251               1.487610               1.566271 
##           Foot.Lam.num           Hand.Lam.num             Avg.lnSVL2 
##               1.616464               1.652438               1.519654 
##              Avg.ln.t1 
##               1.200938

It really is that easy. For Pagel's λ the function will return both a fitted λ value, as well as a log-likelihood. The result is that our apply( ) call produces a list which could require some post-processing. Fortunately, we can use an R base function like simplify2array to get our results in a more compact format:

lambda<-t(simplify2array(apply(Anolis.data,2,phylosig,
    tree=Anolis.tree,method="lambda")))
lambda

##                        lambda   logL     
## AVG.SVL                1.016502 -3.810016
## AVG.hl                 1.018967 -7.919063
## AVG.hw                 1.01674  -18.27221
## AVG.hh                 1.007532 -23.00704
## AVG.ljl                1.016123 -10.03831
## AVG.outlever           1.018473 -9.490815
## AVG.jugal.to.symphysis 1.021371 -6.919674
## AVG.femur              1.013122 -11.72327
## AVG.tibia              1.021411 -7.14409 
## AVG.met                1.021765 -5.985474
## AVG.ltoe.IV            1.021493 -13.6596 
## AVG.toe.IV.lam.width   1.001934 -29.00784
## AVG.humerus            1.011935 -11.56159
## AVG.radius             1.019787 -12.34023
## AVG.lfing.IV           1.00934  -17.31632
## AVG.fing.IV.lam.width  1.008937 -32.1535 
## AVG.pelv.ht            1.017735 -22.90027
## AVG.pelv.wd            1.021306 -9.937602
## Foot.Lam.num           1.020469 55.17366 
## Hand.Lam.num           1.020239 44.33899 
## Avg.lnSVL2             1.016559 -6.030245
## Avg.ln.t1              1.014146 -24.90667

Pretty straightforward, no?

Note that if we had just send a column of the data frame to phylosig( ) as follows:

phylosig(Anolis.tree,Anolis.data[,"AVG.met"])

## [1] "x has no names; assuming x is in the same order as tree$tip.label"

## [1] 0.2514011

that wouldn't do because the order of the rows in our data frame do not correspond to the order of the tip labels in the tree. We can resolve this by using the function setNames:

phylosig(Anolis.tree,setNames(Anolis.data[,"AVG.met"],
    rownames(Anolis.data)))

## [1] 1.585221

but I like our apply( ) calls much better.

logLik method for brownie.lite

I just pushed a small update to phytools giving the popular brownie.lite function a logLik S3 method.

brownie.lite fits a multi-rate Brownian motion model originally developed by Brian O'Meara et al. (2006).

Here, I'll show how it can be used in combination with the generic AIC method to do model selection or to compute model weights:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## The tree includes a mapped, 2-state discrete character with states:
##  a, b
## 
## Rooted; includes branch lengths.

x

##            A            B            C            D            E 
##  -4.42905986  -4.33868707  -5.30046141  -3.05157557  -3.82136084 
##            F            G            H            I            J 
##  -3.84278970  -4.28031068  -4.79455901  -5.81169561  -6.34045086 
##            K            L            M            N            O 
##  -1.65886570  -9.33620130 -10.28889806   0.41905850   0.62457524 
##            P            Q            R            S            T 
##   1.69227697   1.14144175  -0.02283264  -1.98316222  -3.61491689 
##            U            V            W            X            Y 
##   1.71643871   2.75486352   0.13485761   0.21096964   0.47606198 
##            Z 
##  -0.38421995

cols<-setNames(c("blue","red"),c("a","b"))
phenogram(tree,x,colors=cols)

##            1            2            3            4            5 
##  -5.01570904  -4.67036911  -5.70806984  -2.94258054  -3.63500326 
##            6            7            8            9           10 
##  -3.98014777  -4.32528265  -5.36044268  -6.05326441  -6.58088057 
##           11           12           13           14           15 
##  -1.76580067  -9.37186885 -10.28889806   0.36848651   1.05920932 
##           16           17           18           19           20 
##   1.74996436   1.40415384  -0.66609199  -2.11198563  -3.28891866 
##           21           22           23           24           25 
##   2.09536923   2.75486352  -0.32092830   0.02392802   0.71404775 
##           26 
##  -1.01180397

## fit two-rate model using brownie.lite
fit<-brownie.lite(tree,x)
fit

## ML single-rate model:
##  s^2 se  a   k   logL
## value    2.2226  0.6164  -1.4325 2   -45.13  
## 
## ML multi-rate model:
##  s^2(a)  se(a)   s^2(b)  se(b)   a   k   logL    
## value    1.0376  0.4312  3.6835  1.5578  -0.6666 3   -43.2547
## 
## P-value (based on X^2): 0.0528 
## 
## R thinks it has found the ML solution.

liks<-logLik(fit)
liks

## single-rate  multi-rate 
##   -45.12998   -43.25468 
## attr(,"df")
## [1] 2 3

AIC<-setNames(AIC(fit),names(liks))
AIC

## single-rate  multi-rate 
##    94.25996    92.50935

aic.w(AIC)

## single-rate  multi-rate 
##    0.294152    0.705848

That could be helpful.

The data were simulated as follows:

set.seed(1)
tree<-pbtree(n=26,tip.label=LETTERS)
tree<-paintSubTree(tree,29,"b","a",stem=0.5)
x<-sim.rates(tree,setNames(c(1,2.5),letters[1:2]))