Updates to phylosig()

I made a few updates to my function phylosig(), which computes phylogenetic signal using two different methods (method="K" and method="lambda"); and can also conduct hypothesis tests of the null hypothesis that phylogenetic signal is absent. I have posted this function on my R-phylogenetics page (direct link to code here).

Most of the updates are designed to allow the function to accept data vectors with missing data (either as missing elements in the vector, or as NA); as well as data vectors that contain data for species not represented in the tree. In the former case (where we have tips in the tree that are not in the data vector), I use drop.tip() from the {ape} package to remove these tips; and in the latter case (where we have data for tips not present in the tree), we simply excise these data elements from our vector. In both cases, an informative message is triggered. Thanks to Christofer's comments on my previous post for inspiring these improvements.

The other update allows the function to return lambda>1.0. There was some discussion about the theoretical maximum for lambda on the R-sig-phylo email listserve (e.g., here and subsequent messages). I now first compute the maximum possible value for lambda (this is the maximum value for lambda for which the likelihood expression is still defined), as follows:

maxLambda<-max(C)/max(C[upper.tri(C)])

And then I optimize lambda for the interval (0,maxLambda):

res<-optimize(f=likelihood,interval=c(0,maxLambda),y=x,C=C,maximum=TRUE)

This seems to work fine. [Note that I believe that the likelihood expression will be defined from -maxLambda to maxLambda, so if we wanted to allow the possibility of lambda less than 0 we could easily add that here.]