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.
Hi,
ReplyDeleteWhere can I find the likMlambda function?
It is in the current CRAN version of phytools.
Delete