Tuesday, February 11, 2014

Small fix in evol.vcv and new S3 print method

A phytools user pointed out that there was a peculiar 'bug' in the phytools function evol.vcv due to an inconsistency between the variable names in the function definition for an internally used function and the names of the variables used inside this function. It is peculiar because (by a stroke of luck) the internally used variable name is already defined correctly within the main function which means that this bug can never cause a problem. I have fixed it anyway in a new function version and phytools build (here).

Since I already had the source code for this function open, I decided to add an S3 print method for objects returned by evol.vcv. This is in the style of what I'd already done for brownie.lite. This seems to work pretty well. Here's a quick demo.

> require(phytools)
Loading required package: phytools
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.3.91’
> tree<-pbtree(n=26,tip.label=LETTERS[26:1])
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-c("A","B")
> tree<-sim.history(tree,Q,anc="A")
> plotSimmap(tree,lwd=3)
no colors provided. using the following legend:
      A       B
"black"   "red"
> V<-list(matrix(c(1,0,0,1),2,2),matrix(c(1,1.2,1.2,2),2,2))
> names(V)<-c("A","B")
> V
$A
     [,1] [,2]
[1,]    1    0
[2,]    0    1

$B
     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  2.0

> X<-sim.corrs(tree,vcv=V)
> fit<-evol.vcv(tree,X)
> fit
ML single-matrix model:
        R[1,1]  R[1,2]  R[2,2]  k       log(L)
fitted  0.675   0.355   1.5944  5       -68.1635
      
ML multi-matrix model:
        R[1,1]  R[1,2]  R[2,2]  k       log(L)
A       0.7789  -0.3202 0.4378  8       -64.6364
B       0.6622  1.0262  2.7149

P-value (based on X^2): 0.0702

R thinks it has found the ML solution.

This is pretty much what we were going for. It might get a little messy if we have a lot of columns in X. That's it.

No comments:

Post a Comment