anova for "fitPagel" object class in phytools

I just added a S3 generic anova method to phytools for the object class "fitPagel".

This is designed to facilitate comparing results across different variants of the Pagel (1994) correlated binary trait evolution model – such as, for instance, when we model that trait x dependents on trait y, but not the converse.

Here's a quick demo using a dataset for paternal care and spawning mode in bony fishes from Benun Sutton & Wilson (2019).

First, read in and plot our tree & data.

library(phytools)
packageVersion("phytools")

## [1] '1.1.13'

bonyfish.tree<-read.tree(file="http://www.phytools.org/Rbook/7/bonyfish.tre")
bonyfish.data<-read.csv(file="http://www.phytools.org/Rbook/7/bonyfish.csv",
    row.names=1,stringsAsFactors=TRUE)
object<-plotTree.datamatrix(bonyfish.tree,bonyfish.data,
    fsize=0.5,yexp=1,header=FALSE,xexp=1.45,
    palettes=c("YlOrRd","PuBuGn"))
leg<-legend(x="topright",names(object$colors$spawning_mode),
    cex=0.7,pch=22,pt.bg=object$colors$spawning_mode,
    pt.cex=1.5,bty="n",title="spawning mode")
leg<-legend(x=leg$rect$left+4.7,y=leg$rect$top-leg$rect$h,
    names(object$colors$paternal_care),cex=0.7,
    pch=22,pt.bg=object$colors$paternal_care,pt.cex=1.5,
    bty="n",title="paternal care")

Next, pull out the two binary characters of interest.

spawning_mode<-setNames(bonyfish.data[,1],
    rownames(bonyfish.data))
head(spawning_mode)

##         Xenomystus_nigri       Chirocentrus_dorab     Talismania_bifurcata Alepocephalus_tenebrosus 
##                     pair                    group                    group                    group 
##     Misgurnus_bipartitus     Opsariichthys_bidens 
##                     pair                     pair 
## Levels: group pair

paternal_care<-setNames(bonyfish.data[,2],
    rownames(bonyfish.data))
head(paternal_care)

##         Xenomystus_nigri       Chirocentrus_dorab     Talismania_bifurcata Alepocephalus_tenebrosus 
##                     male                     none                     none                     none 
##     Misgurnus_bipartitus     Opsariichthys_bidens 
##                     none                     none 
## Levels: male none

Next, I'm going to run fitPagel three times.

First, I'll fit a model in which x (spawning mode) depends on y (paternal care), but not the reverse.

Second, I'll fit a model in which paternal care depends on spawning mode.

Finally, I'll fit a model in which both x depends on y, and y depends on x.

Remember, with Pagel's (1994) model when we say “dependent” – what we mean is that the evolutionary rates of change depend for one trait are affected by the state of the other.

dep.spawning_mode<-fitPagel(bonyfish.tree,x=spawning_mode,
    y=paternal_care,dep.var="x",pi="fitzjohn")
dep.paternal_care<-fitPagel(bonyfish.tree,x=spawning_mode,
    y=paternal_care,dep.var="y",pi="fitzjohn")
dep.xy<-fitPagel(bonyfish.tree,x=spawning_mode,
    y=paternal_care,dep.var="xy",pi="fitzjohn")

Now, we can use the anova method as follows.

anova(dep.spawning_mode,dep.paternal_care,
    dep.xy)

##                      log(L) d.f.      AIC     weight
## independent       -60.78126    4 129.5625 0.02977125
## dep.spawning_mode -55.64916    6 123.2983 0.68241595
## dep.paternal_care -57.53577    6 127.0715 0.10344413
## dep.xy            -54.95787    8 125.9157 0.18436867

This shows us that the best-supported model is one in which spawning mode depends on paternal care, but not the reverse.

Let's see the plot of this model.

plot(dep.paternal_care,lwd.by.rate=TRUE)

Pretty neat.

Note that anova also works for just a single fitPagel model object, in which case it merely extracts the log-likelihoods, AIC, and Akaike weights of our independent & dependent models.

For instance.

anova(dep.xy)

##                log(L) d.f.      AIC    weight
## independent -60.78126    4 129.5625 0.1390271
## dep.xy      -54.95787    8 125.9157 0.8609729

That's all folks.