Today I wrote an S3 print method for objects returned by the phytools function brownie.lite. This function implements the method of O'Meara et al. (2006) in which we fit different rates of evolution to different parts of a phylogeny (sometimes, say, determined by a mapped discrete character). Print methods are nice because they allow us to tell R how to summarize the content of a complicated object, like a phylogeny or the results of numerical optimization.
The object returned by brownie.lite is unchanged, except that it is now assigned the class attribute "brownie.lite". Here's how it works:
> require(phytools)
> packageVersion("phytools")
[1] ‘0.3.83’
> ## first let's simulate some data
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-letters[1:2]
> tree<-sim.history(pbtree(n=100,scale=1),Q,anc="a")
> ## simulate with a rate difference among states
> x<-sim.rates(tree,setNames(c(1,2),letters[1:2]))
> ## simulate without a rate difference
> y<-fastBM(tree)
> plotSimmap(tree,setNames(c("blue","red"),letters[1:2]), lwd=3,ftype="off")
> packageVersion("phytools")
[1] ‘0.3.83’
> ## first let's simulate some data
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-letters[1:2]
> tree<-sim.history(pbtree(n=100,scale=1),Q,anc="a")
> ## simulate with a rate difference among states
> x<-sim.rates(tree,setNames(c(1,2),letters[1:2]))
> ## simulate without a rate difference
> y<-fastBM(tree)
> plotSimmap(tree,setNames(c("blue","red"),letters[1:2]), lwd=3,ftype="off")
> fitx<-brownie.lite(tree,x)
> fitx # this is the same as print(fitx)
ML single-rate model:
s^2 se a k logL
parameter 1.5484 0.2191 -0.4682 2 -73.73
ML multi-rate model:
s^2(a) se(a) s^2(b) se(b) a k logL
parameter 0.732 0.2243 1.8163 0.2955 -0.499 3 -70.93
P-value (based on X^2): 0.0179
R thinks it has found the ML solution.
> fity<-brownie.lite(tree,y)
> fity
ML single-rate model:
s^2 se a k logL
parameter 1.0646 0.1507 -0.7218 2 -55.00
ML multi-rate model:
s^2(a) se(a) s^2(b) se(b) a k logL
parameter 0.8925 0.2661 1.1209 0.1841 -0.7271 3 -54.79
P-value (based on X^2): 0.5184
R thinks it has found the ML solution.
> fitx # this is the same as print(fitx)
ML single-rate model:
s^2 se a k logL
parameter 1.5484 0.2191 -0.4682 2 -73.73
ML multi-rate model:
s^2(a) se(a) s^2(b) se(b) a k logL
parameter 0.732 0.2243 1.8163 0.2955 -0.499 3 -70.93
P-value (based on X^2): 0.0179
R thinks it has found the ML solution.
> fity<-brownie.lite(tree,y)
> fity
ML single-rate model:
s^2 se a k logL
parameter 1.0646 0.1507 -0.7218 2 -55.00
ML multi-rate model:
s^2(a) se(a) s^2(b) se(b) a k logL
parameter 0.8925 0.2661 1.1209 0.1841 -0.7271 3 -54.79
P-value (based on X^2): 0.5184
R thinks it has found the ML solution.
This version of phytools can be downloaded here and installed from source.
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.