Sunday, May 29, 2016

Some cool additional features for "fitPagel" plot method

I just added a whole suite of additional features to the S3 plotting method for objects of class "fitPagel" that I first described earlier today.

Among these, perhaps the most interesting is the addition of a logical argument lwd.by.rate which allows the user to plot the transition matrices in which the arrow thickness is proportional to the rate. I also added features to control the headings (main, which can be a vector of two elements), heading font size (cex.main), subheading font size (cex.sub), trait label font size (cex.traits), and rate font size (cex.rates). In addition, I added a feature to split the traits across two lines in the dependent model should the summed number of character for the two traits be larger than, if I remember correctly, five.

Here's a demo. The data are simulated, but the consist of two binary trait vectors, x & y. I'm going to fit three models: one in which x depends on y, a second that is the converse, and a third in which x and y depend each on the other. After each model, I will plot the results.

library(phytools)
tree
## 
## Phylogenetic tree with 300 tips and 299 internal nodes.
## 
## Tip labels:
##  t149, t150, t164, t165, t59, t60, ...
## 
## Rooted; includes branch lengths.
head(x,n=20)
##          t1          t2          t3          t4          t5          t6 
##     aquatic terrestrial     aquatic terrestrial     aquatic terrestrial 
##          t7          t8          t9         t10         t11         t12 
## terrestrial     aquatic terrestrial     aquatic     aquatic terrestrial 
##         t13         t14         t15         t16         t17         t18 
## terrestrial     aquatic terrestrial terrestrial terrestrial terrestrial 
##         t19         t20 
## terrestrial terrestrial 
## Levels: terrestrial aquatic
head(y,n=20)
##     t1     t2     t3     t4     t5     t6     t7     t8     t9    t10 
## keeled   flat keeled   flat keeled   flat   flat keeled   flat keeled 
##    t11    t12    t13    t14    t15    t16    t17    t18    t19    t20 
## keeled   flat   flat keeled   flat keeled   flat   flat   flat   flat 
## Levels: flat keeled
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:
##                    terrestrial|flat terrestrial|keeled aquatic|flat
## terrestrial|flat         -1.6338807          0.7929193    0.8409613
## terrestrial|keeled        0.4595752         -1.3005366    0.0000000
## aquatic|flat              0.5940199          0.0000000   -1.3869392
## aquatic|keeled            0.0000000          0.5940199    0.4595752
##                    aquatic|keeled
## terrestrial|flat        0.0000000
## terrestrial|keeled      0.8409613
## aquatic|flat            0.7929193
## aquatic|keeled         -1.0535951
## 
## Dependent (x only) model rate matrix:
##                    terrestrial|flat terrestrial|keeled aquatic|flat
## terrestrial|flat         -1.1647920          0.7981942    0.3665978
## terrestrial|keeled        0.4506319         -5.1481648    0.0000000
## aquatic|flat              3.0144733          0.0000000   -3.8126675
## aquatic|keeled            0.0000000          0.3205245    0.4506319
##                    aquatic|keeled
## terrestrial|flat        0.0000000
## terrestrial|keeled      4.6975329
## aquatic|flat            0.7981942
## aquatic|keeled         -0.7711564
## 
## 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,
    main=c("a) Independent model",
    "b) Dependent (terrestrial | aquatic) model"))

plot of chunk unnamed-chunk-1

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:
##                    terrestrial|flat terrestrial|keeled aquatic|flat
## terrestrial|flat         -1.6338807          0.7929193    0.8409613
## terrestrial|keeled        0.4595752         -1.3005366    0.0000000
## aquatic|flat              0.5940199          0.0000000   -1.3869392
## aquatic|keeled            0.0000000          0.5940199    0.4595752
##                    aquatic|keeled
## terrestrial|flat        0.0000000
## terrestrial|keeled      0.8409613
## aquatic|flat            0.7929193
## aquatic|keeled         -1.0535951
## 
## Dependent (y only) model rate matrix:
##                    terrestrial|flat terrestrial|keeled aquatic|flat
## terrestrial|flat         -1.2400499          0.3923515    0.8476983
## terrestrial|keeled        3.1664689         -4.0141672    0.0000000
## aquatic|flat              0.6916232          0.0000000   -4.6332749
## aquatic|keeled            0.0000000          0.6916232    0.1515701
##                    aquatic|keeled
## terrestrial|flat        0.0000000
## terrestrial|keeled      0.8476983
## aquatic|flat            3.9416517
## aquatic|keeled         -0.8431933
## 
## 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,
    main=c("a) Independent model",
    "b) Dependent (flat | keeled) model"))

plot of chunk unnamed-chunk-1

fit.xy<-fitPagel(tree,x,y,dep.var="xy")
fit.xy
## 
## Pagel's binary character correlation test:
## 
## Assumes "ARD" substitution model for both characters
## 
## Independent model rate matrix:
##                    terrestrial|flat terrestrial|keeled aquatic|flat
## terrestrial|flat         -1.6338807          0.7929193    0.8409613
## terrestrial|keeled        0.4595752         -1.3005366    0.0000000
## aquatic|flat              0.5940199          0.0000000   -1.3869392
## aquatic|keeled            0.0000000          0.5940199    0.4595752
##                    aquatic|keeled
## terrestrial|flat        0.0000000
## terrestrial|keeled      0.8409613
## aquatic|flat            0.7929193
## aquatic|keeled         -1.0535951
## 
## Dependent (x & y) model rate matrix:
##                    terrestrial|flat terrestrial|keeled aquatic|flat
## terrestrial|flat         -1.0668076          0.6538038    0.4130038
## terrestrial|keeled        0.9094271         -4.6819194    0.0000000
## aquatic|flat              2.6064761          0.0000000   -4.2118594
## aquatic|keeled            0.0000000          0.3624449    0.3440521
##                    aquatic|keeled
## terrestrial|flat        0.0000000
## terrestrial|keeled      3.7724923
## aquatic|flat            1.6053833
## aquatic|keeled         -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(fit.xy,lwd.by.rate=TRUE,
    main=c("a) Independent model",
    "b) Full dependent model"))

plot of chunk unnamed-chunk-1

Pretty neat.

1 comment:

  1. Is there a way to use lwd.by.rate=TRUE for plotting fitMk outputs?? I can't figure it out, or find any examples.
    Thanks!

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.