S3 plot method for objects of class "fitPagel"

I just added an S3 plot method for objects of class "fitPagel".

A "fitPagel" object normally results from a call to the function "fitPagel" which (sensibly) fits Pagel's (1994) model to test whether the evolution of two binary characters is correlated. (That is, if the state of one character influences the rate of a second, or vice versa.)

Probably the best way to explain the plot method is to demonstrate it.

For that I will use data simulated for another exercise. Here I can read this data in from the web as I have used them in various workshops:

library(phytools)

## Loading required package: ape

## Loading required package: maps

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

packageVersion("phytools")

## [1] '0.5.33'

tree<-read.tree("http://www.phytools.org/Barcelona2016/data/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("http://www.phytools.org/Barcelona2016/data/fitPagel-data.csv",
    row.names=1)
head(X) ## so you can see the data format

##    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

## pull out first two columns as factors:
x<-setNames(X[,1],rownames(X))
y<-setNames(X[,2],rownames(X))
obj<-fitPagel(tree,x,y)
obj

## 
## 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.6338807  0.7929193  0.8409613  0.0000000
## a|b  0.4595752 -1.3005366  0.0000000  0.8409613
## b|a  0.5940199  0.0000000 -1.3869392  0.7929193
## b|b  0.0000000  0.5940199  0.4595752 -1.0535951
## 
## Dependent (x & y) model rate matrix:
##            a|a        a|b        b|a        b|b
## a|a -1.0668076  0.6538038  0.4130038  0.0000000
## a|b  0.9094271 -4.6819194  0.0000000  3.7724923
## b|a  2.6064761  0.0000000 -4.2118594  1.6053833
## b|b  0.0000000  0.3624449  0.3440521 -0.7064971
## 
## 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

plot(obj)

This is rather crude right now - it doesn't adjust the spacing of the arrows to accommodate different state lengths. Nonetheless, I think it looks cool. Please let me know if you encounter any bugs.

The plot method code is as follows:

plot.fitPagel<-function(x,...){
    if(hasArg(signif)) signif<-list(...)$signif
    else signif<-3
    par(mfrow=c(2,1))
    ## INDEPENDENT MODEL
    plot.new()
    par(mar=c(1.1,2.1,3.1,2.1))
    plot.window(xlim=c(0,2),ylim=c(0,1),asp=1)
    mtext("a) Independent model",side=3,adj=0,line=1.2,cex=1.2)
    ## trait 1
    text(x=0.15,y=1,"Trait 1:")
    arrows(x0=0.4,y0=0.15,y1=0.85,lwd=2,length=0.1)
    arrows(x0=0.45,y0=0.85,y1=0.15,lwd=2,length=0.1)
    text(x=0.425,y=0.95,strsplit(rownames(x$dependent.Q)[1],"|")[[1]][1])
    text(x=0.425,y=0.05,strsplit(rownames(x$dependent.Q)[3],"|")[[1]][1])
    text(x=0.50,y=0.5,round(x$independent.Q[1,3],signif),cex=0.8,srt=90)
    text(x=0.35,y=0.5,round(x$independent.Q[3,1],signif),cex=0.8,srt=90)
    ## trait 2
    text(x=1.3,y=1,"Trait 2:")  
    arrows(x0=1.55,y0=0.15,y1=0.85,lwd=2,length=0.1)
    arrows(x0=1.60,y0=0.85,y1=0.15,lwd=2,length=0.1)
    text(x=1.575,y=0.95,strsplit(rownames(x$dependent.Q)[1],"|")[[1]][3])
    text(x=1.575,y=0.05,strsplit(rownames(x$dependent.Q)[2],"|")[[1]][3])
    text(x=1.65,y=0.5,round(x$independent.Q[1,2],signif),cex=0.8,srt=90)
    text(x=1.50,y=0.5,round(x$independent.Q[2,1],signif),cex=0.8,srt=90)
    ## DEPENDENT MODEL
    plot.new()
    par(mar=c(1.1,2.1,3.1,2.1))
    plot.window(xlim=c(0,2),ylim=c(0,1),asp=1)
    mtext("b) Dependent model",side=3,adj=0,line=1.2,cex=1.2)
    arrows(x0=0.15,y0=0.15,y1=0.85,lwd=2,length=0.1)
    arrows(x0=0.2,y0=0.85,y1=0.15,lwd=2,length=0.1)
    arrows(x0=1.6,y0=0.05,x1=0.4,lwd=2,length=0.1)
    arrows(x0=0.4,y0=0.1,x1=1.6,lwd=2,length=0.1)
    arrows(x0=1.8,y0=0.15,y1=0.85,lwd=2,length=0.1)
    arrows(x0=1.85,y0=0.85,y1=0.15,lwd=2,length=0.1)
    arrows(x0=1.6,y0=0.9,x1=0.4,lwd=2,length=0.1)
    arrows(x0=0.4,y0=0.95,x1=1.6,lwd=2,length=0.1)
    ## add states
    text(x=0.175,y=0.95,
        paste(strsplit(rownames(x$dependent.Q)[1],"|")[[1]][c(1,3)],
        collapse=", "))
    text(x=1.825,y=0.95,
        paste(strsplit(rownames(x$dependent.Q)[2],"|")[[1]][c(1,3)],
        collapse=", "))
    text(x=1.825,y=0.05,
        paste(strsplit(rownames(x$dependent.Q)[4],"|")[[1]][c(1,3)],
        collapse=", "))
    text(x=0.175,y=0.05,
        paste(strsplit(rownames(x$dependent.Q)[3],"|")[[1]][c(1,3)],
        collapse=", "))
    ## add rates
    text(x=1,y=1,round(x$dependent.Q[1,2],signif),cex=0.8)
    text(x=1,y=0.85,round(x$dependent.Q[2,1],signif),cex=0.8)
    text(x=1.9,y=0.5,round(x$dependent.Q[2,4],signif),cex=0.8,srt=90)
    text(x=1.75,y=0.5,round(x$dependent.Q[4,2],signif),cex=0.8,srt=90)
    text(x=1,y=0,round(x$dependent.Q[4,3],signif),cex=0.8)
    text(x=1,y=0.15,round(x$dependent.Q[3,4],signif),cex=0.8)
    text(x=0.1,y=0.5,round(x$dependent.Q[3,1],signif),cex=0.8,srt=90)
    text(x=0.25,y=0.5,round(x$dependent.Q[1,3],signif),cex=0.8,srt=90)
}

The latest version of phytools can be installed from GitHub using devtools:

library(devtools)
install_github("liamrevell/phytools",dependencies=FALSE)

That's it.