Tuesday, September 11, 2012

Joint estimation of Pagel's λ for multiple characters

In a previous post, I described adding a previously internal function for computing the likelihood of a single value of λ for multiple traits to the NAMESPACE of phytools. This function (likMlambda) can now be used in a very simple way to find the ML joint λ estimate for multiple characters. This is described (if I remember correctly) in Freckleton et al. (2002) and implemented internally in the phytools functions phyl.cca & phyl.pca, but not in any stand alone functions. Here is the function:

joint.lambda<-function(tree,X){
  result<-optimize(f=likMlambda,X=X,C=vcv(tree),interval=c(0,1),     maximum=TRUE)
  return(list(lambda=result$maximum,logL=result$objective[1]))
}


Let's try it out:

> tree<-pbtree(n=100)
> X<-fastBM(lambdaTree(tree,0.7),nsim=10)
> joint.lambda(tree,X)
$lambda
[1] 0.6722762
$logL
[1] -1816.575


Pretty cool.

2 comments:

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.