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"))
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"))
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"))
Pretty neat.
Is there a way to use lwd.by.rate=TRUE for plotting fitMk outputs?? I can't figure it out, or find any examples.
ReplyDeleteThanks!