Matt Johnson, a graduate student at Duke, just identified an interesting problem with the phytools function phylosig. Basically, he was finding that phylosig(...,method="lambda",se=se) was giving discrepant results when he rescaled the branch lengths of the tree to be very short (say, scaling all edge lengths by 1/1000). This was a problem unique to the optimization of λ with sampling error in the estimation of species means. This immediately pointed towards the R base package multidimensional optimizer, optim, because no other methods in phylosig use optim aside from λ with sampling error.
To reproduce the error that he found, Matt suggested something along the lines of the following code:
> require(phytools); require(geiger)
> se<-sqrt(rchisq(n=200,df=2)); names(se)<-tree$tip.label
Ok, so far so good. Now let's try and rescale the tree to have really short branches. In theory this should not affect our estimate of phylogenetic signal.
> tree1000<-tree; tree1000$edge.length<-tree$edge.length/1000
Obviously, there's some kind of problem here. What I've been able to figure out is that it is a convergence issue - for some reason, convergence is more difficult for the larger value of σ2 that results from rescaling the branch lengths to have very short lengths. [Note that for any change to the branch lengths of the tree by factor 1/k, the fitted evolutionary rate should increase by factor k.] I've made some headway towards resolving this issue and will post more, along with an updated version of phylosig shortly.
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.