R function for Pagel's 1994 correlation method

I just posted some code to fit the Pagel (1994) method for detecting “correlated” evolution of two binary traits. I say correlated because it might be more appropriate to call this a dependence - i.e., the rate of character 2 depends on 1 & vice versa.

This is actually quite simple to do because the test is merely a test of, for a two character combination written as [01] (meaning state 0 for the first character & state 1 for the second), the rates of transition from [00]->[01] equal [10]->[11]; [00]->[10] equal [01]->[11]; etc. For this we can just paste our two input characters together, fit both a model in which the dependence exists and one in which it does not, and then compare them. The dependence model has twice as many transition rates as the independence model, so this model will have (for two binary traits) 4 more parameters.

Note that I accomplished this by borrowing ace in the ape package or fitDiscrete in the geiger package to actually fit the models. All my code really does is pre- & post-process the data & results respectively!

Here's a demo:

## first load packages & source code
library(phytools)
library(geiger)
source("fitPagel.R")
.check.pkg<-phytools:::.check.pkg
## now let's simulate some uncorrelated data
tree<-pbtree(n=300,scale=1)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
tt1<-sim.history(tree,Q)

## Done simulation(s).

tt2<-sim.history(tree,Q)

## Done simulation(s).

## these are uncorrelated, see:
par(mfrow=c(1,2))
plotSimmap(tt1,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1)
plotSimmap(tt2,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1,direction="leftwards")

x<-tt1$states
y<-tt2$states
fit.ape<-fitPagel(tree,x,y)
fit.ape

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.9712  1.0896  0.8816  0.0000
## a|b  0.9774 -1.8589  0.0000  0.8816
## b|a  0.7531  0.0000 -1.8426  1.0896
## b|b  0.0000  0.7531  0.9774 -1.7304
## 
## Dependent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.5606  0.8578  0.7028  0.0000
## a|b  0.9261 -1.8634  0.0000  0.9373
## b|a  0.3652  0.0000 -1.5764  1.2112
## b|b  0.0000  1.2754  1.0567 -2.3321
## 
## Model fit:
##             log-likelihood
## independent         -238.7
## dependent           -237.5
## 
## Hypothesis test result:
##   likelihood-ratio:  2.390211 
##   p-value:  0.6643971

fit.geiger<-fitPagel(tree,x,y,method="fitDiscrete")

## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41
## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41

fit.geiger

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.8847  0.9942  0.8906  0.0000
## a|b  0.9699 -1.8605  0.0000  0.8906
## b|a  0.7112  0.0000 -1.7053  0.9942
## b|b  0.0000  0.7112  0.9699 -1.6811
## 
## Dependent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.5000  0.8192  0.6809  0.0000
## a|b  0.8987 -1.8596  0.0000  0.9609
## b|a  0.3704  0.0000 -1.4820  1.1116
## b|b  0.0000  1.1176  1.1075 -2.2251
## 
## Model fit:
##             log-likelihood
## independent         -239.3
## dependent           -238.3
## 
## Hypothesis test result:
##   likelihood-ratio:  2.096727 
##   p-value:  0.7179738

OK, now let's try correlated data. To do that I will have to build a big transition matrix myself for simulation!

Q<-matrix(c(0,0.5,0.5,0,2,0,0,2,2,0,0,2,0,0.5,0.5,0),4,4,byrow=TRUE)
rownames(Q)<-colnames(Q)<-c("aa","ab","ba","bb")
diag(Q)<--rowSums(Q)
tt<-sim.history(tree,t(Q))

## Note - the rate of substitution from i->j should be given by Q[j,i].
## Done simulation(s).

tt1<-mergeMappedStates(tt,c("aa","ab"),"a")
tt1<-mergeMappedStates(tt1,c("ba","bb"),"b")
tt2<-mergeMappedStates(tt,c("aa","ba"),"a")
tt2<-mergeMappedStates(tt2,c("ab","bb"),"b")
## these data are correlated, see:
par(mfrow=c(1,2))
plotSimmap(tt1,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1)
plotSimmap(tt2,setNames(c("blue","red"),letters[1:2]),ftype="off",lwd=1,direction="leftwards")

x<-getStates(tt1,"tips")
y<-getStates(tt2,"tips")
fit.ape<-fitPagel(tree,x,y)
fit.ape

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a     a|b     b|a     b|b
## a|a -1.4325  0.9831  0.4493  0.0000
## a|b  0.9593 -1.4086  0.0000  0.4493
## b|a  1.0401  0.0000 -2.0233  0.9831
## b|b  0.0000  1.0401  0.9593 -1.9994
## 
## Dependent model rate matrix:
##        a|a     a|b     b|a    b|b
## a|a -1.043  0.8539  0.1893  0.000
## a|b  5.543 -7.9019  0.0000  2.359
## b|a  1.202  0.0000 -4.4114  3.209
## b|b  0.000  1.0849  0.4054 -1.490
## 
## Model fit:
##             log-likelihood
## independent         -241.1
## dependent           -207.8
## 
## Hypothesis test result:
##   likelihood-ratio:  66.63346 
##   p-value:  1.164753e-13

fit.geiger<-fitPagel(tree,x,y,method="fitDiscrete")

## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41
## Warning: Parameter estimates appear at bounds:
##  q14
##  q23
##  q32
##  q41

fit.geiger

## 
##   Pagel's binary character correlation test:
## 
## Indepedent model rate matrix:
##         a|a    a|b     b|a     b|b
## a|a -1.3892  1.006  0.3832  0.0000
## a|b  0.9236 -1.307  0.0000  0.3832
## b|a  1.0639  0.000 -2.0699  1.0061
## b|b  0.0000  1.064  0.9236 -1.9875
## 
## Dependent model rate matrix:
##         a|a     a|b        b|a    b|b
## a|a -0.8951  0.8951  1.799e-16  0.000
## a|b  5.4502 -7.7907  0.000e+00  2.340
## b|a  1.9288  0.0000 -3.958e+00  2.029
## b|b  0.0000  1.0086  4.105e-01 -1.419
## 
## Model fit:
##             log-likelihood
## independent         -241.9
## dependent           -208.8
## 
## Hypothesis test result:
##   likelihood-ratio:  66.16744 
##   p-value:  1.460391e-13

That's it, I guess. This is brand new so I have no idea how it works but feedback is welcome. I'm also well aware of the recent critique of this method by Maddision & Fitzjohn.