Exercise 9: Testing for an evolutionary correlation between binary characters using Pagel’s (1994) method

This tutorial is a short exercise designed to fit and test Pagel's (1994) method for testing for an evolutionary relationship between two binary characters.

To do this exercise, you will first need to download the following dataset & tree, although they can also be loaded directly from the web:

  1. fitPagel-data.csv
  2. fitPagel-tree.tre.

Let's load phytools and then our data and tree:

library(phytools)
tree<-read.tree("fitPagel-tree.tre")
tree
## 
## Phylogenetic tree with 300 tips and 299 internal nodes.
## 
## Tip labels:
##  t149, t150, t164, t165, t59, t60, ...
## 
## Rooted; includes branch lengths.
X<-read.csv("fitPagel-data.csv",row.names=1,header=TRUE)
## X is a big matrix, so let's just visualize the top part
dim(X)
## [1] 300   3
head(X)
##    x y z
## t1 b b a
## t2 a a a
## t3 b b b
## t4 a a a
## t5 b b a
## t6 a a a

Next, let's pull out the characters we're going to analyze as vectors:

x<-setNames(X$x,rownames(X))
y<-setNames(X$y,rownames(X))
z<-setNames(X$z,rownames(X))

OK, having done this, we will use the phytools function fitPagel to fit Pagel's (1994) model. To get a sense of our data before we do, we should first visualize the data on the tree. Let's try this:

par(mfrow=c(1,2))
plot(tree,show.tip.label=FALSE,no.margin=TRUE)
par(fg="transparent")
tiplabels(pie=to.matrix(x[tree$tip.label],c("a","b")),piecol=c("blue","red"),
    cex=0.3)
par(fg="black")
add.simmap.legend(colors=setNames(c("blue","red"),c("a","b")),prompt=FALSE,
    x=0,y=10,fsize=1)
par(fg="transparent")
plot(tree,show.tip.label=FALSE,no.margin=TRUE,direction="leftwards")
tiplabels(pie=to.matrix(y[tree$tip.label],c("a","b")),piecol=c("blue","red"),
    cex=0.3)

plot of chunk unnamed-chunk-3

par(fg="black")

There may be a pattern of association, but from a visual inspection it's hard to tell.

OK, so we can repeat that also with the third column of X:

par(mfrow=c(1,2))
plot(tree,show.tip.label=FALSE,no.margin=TRUE)
par(fg="transparent")
tiplabels(pie=to.matrix(x[tree$tip.label],c("a","b")),piecol=c("blue","red"),
    cex=0.3)
par(fg="black")
add.simmap.legend(colors=setNames(c("blue","red"),c("a","b")),prompt=FALSE,
    x=0,y=10,fsize=1)
par(fg="transparent")
plot(tree,show.tip.label=FALSE,no.margin=TRUE,direction="leftwards")
tiplabels(pie=to.matrix(z[tree$tip.label],c("a","b")),piecol=c("blue","red"),
    cex=0.3)

plot of chunk unnamed-chunk-4

par(fg="black")

Once again, it's hard to tell whether or not the two characters have evolved in a coordinated fashion.

Now we're ready to run our models. First we will fit a model in which x depends on y, and vice versa; then we will fit a model in which x depends on z and the reverse.

fit.xy<-fitPagel(tree,x,y)
fit.xy
## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.6338771  0.7929164  0.8409607  0.0000000
## a|b  0.4595897 -1.3005505  0.0000000  0.8409607
## b|a  0.5940171  0.0000000 -1.3869335  0.7929164
## b|b  0.0000000  0.5940171  0.4595897 -1.0536068
## 
## Dependent (x & y) model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.0668064  0.6538025  0.4130039  0.0000000
## a|b  0.9094159 -4.6818949  0.0000000  3.7724790
## b|a  2.6064963  0.0000000 -4.2118758  1.6053795
## b|b  0.0000000  0.3624417  0.3440560 -0.7064977
## 
## Model fit:
##             log-likelihood      AIC
## independent      -233.0357 474.0715
## dependent        -206.5951 429.1903
## 
## Hypothesis test result:
##   likelihood-ratio:  52.88121 
##   p-value:  9.023704e-11 
## 
## Model fitting method used was fitMk
fit.xz<-fitPagel(tree,x,z)
fit.xz
## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -2.1937429  1.3527818  0.8409611  0.0000000
## a|b  0.8686789 -1.7096399  0.0000000  0.8409611
## b|a  0.5940193  0.0000000 -1.9468011  1.3527818
## b|b  0.0000000  0.5940193  0.8686789 -1.4626982
## 
## Dependent (x & y) model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -2.1476611  1.7078884  0.4397727  0.0000000
## a|b  0.9879923 -1.9819688  0.0000000  0.9939766
## b|a  0.3351617  0.0000000 -1.0123439  0.6771822
## b|b  0.0000000  0.6617982  0.6723597 -1.3341580
## 
## Model fit:
##             log-likelihood      AIC
## independent      -248.4293 504.8585
## dependent        -246.8280 509.6559
## 
## Hypothesis test result:
##   likelihood-ratio:  3.202592 
##   p-value:  0.5245125 
## 
## Model fitting method used was fitMk

Our results suggest that, indeed, x and y (but not x and z) may be co-evolving. Let's use the S3 plot method to visualize the result:

plot(fit.xy)

plot of chunk unnamed-chunk-6

## or
plot(fit.xy,lwd.by.rate=TRUE)

plot of chunk unnamed-chunk-6

Finally, we can also fit a model in which x depends on y, but not the converse, or y depends on x, but not the converse. We'll focus on x & y only because x and z do not appear to be correlated at all.

To do this we can run the following:

fit.x<-fitPagel(tree,x,y,dep.var="x")
fit.x
## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.6338771  0.7929164  0.8409607  0.0000000
## a|b  0.4595897 -1.3005505  0.0000000  0.8409607
## b|a  0.5940171  0.0000000 -1.3869335  0.7929164
## b|b  0.0000000  0.5940171  0.4595897 -1.0536068
## 
## Dependent (x only) model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.1647923  0.7981948  0.3665975  0.0000000
## a|b  0.4506321 -5.1481654  0.0000000  4.6975333
## b|a  3.0144648  0.0000000 -3.8126596  0.7981948
## b|b  0.0000000  0.3205196  0.4506321 -0.7711517
## 
## Model fit:
##             log-likelihood      AIC
## independent      -233.0357 474.0715
## dependent        -207.6239 427.2478
## 
## Hypothesis test result:
##   likelihood-ratio:  50.8237 
##   p-value:  9.19968e-12 
## 
## Model fitting method used was fitMk
plot(fit.x,lwd.by.rate=TRUE)

plot of chunk unnamed-chunk-7

fit.y<-fitPagel(tree,x,y,dep.var="y")
fit.y
## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.6338771  0.7929164  0.8409607  0.0000000
## a|b  0.4595897 -1.3005505  0.0000000  0.8409607
## b|a  0.5940171  0.0000000 -1.3869335  0.7929164
## b|b  0.0000000  0.5940171  0.4595897 -1.0536068
## 
## Dependent (y only) model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.2400484  0.3923514  0.8476970  0.0000000
## a|b  3.1664563 -4.0141533  0.0000000  0.8476970
## b|a  0.6916106  0.0000000 -4.6332636  3.9416530
## b|b  0.0000000  0.6916106  0.1515679 -0.8431785
## 
## Model fit:
##             log-likelihood      AIC
## independent      -233.0357 474.0715
## dependent        -213.2927 438.5854
## 
## Hypothesis test result:
##   likelihood-ratio:  39.48603 
##   p-value:  2.665125e-09 
## 
## Model fitting method used was fitMk
plot(fit.y,lwd.by.rate=TRUE)

plot of chunk unnamed-chunk-7

Now we have fit four different models for the (co)evolution of x and y. One in which they evolve independently; one in which x depends on y, but not the reverse; one in which y depends on x, but not the reverse; and, finally, one in which x and y are interdependent. We can compare these different possibilities using AIC:

aic<-setNames(c(fit.xy$independent.AIC,
    fit.x$dependent.AIC,
    fit.y$dependent.AIC,
    fit.xy$dependent.AIC),
    c("independent","dependent x",
    "dependent y","dependent x&y"))
aic
##   independent   dependent x   dependent y dependent x&y 
##      474.0715      427.2478      438.5854      429.1903
aic.w(aic)
##   independent   dependent x   dependent y dependent x&y 
##    0.00000000    0.72355637    0.00249763    0.27394600

This shows that the strongest weight of evidence supports the dependent x model, although interdependent x & y is also reasonably strongly supported.

Challenge Problem 6: Detecting the evolutionary correlation between binary traits

Download the following two files:

  1. fitPagel-challenge.tre
  2. fitPagel-challenge.csv

Use the phytools function fitPagel to fit Pagel's model for correlated binary-trait evolution on the tree. What do you find? Visualize the character data on the tree. What does this tell you about the basis for your result?

Written by Liam J. Revell. Last updated 2 August 2017.