In a follow up to the R-sig-phylo thread that initiated yesterday's thread, an R phylogenetics user report that the method fails for high values of α with the error report that the among-species VCV matrix was singularity or nearly so. My guess is that this is due to the optimizer trying values of σ2 which cause the matrix to be numerically indistinguishable from singular. Singular matrices, of course, cannot be inverted, which will cause the likelihood calculations to fail.
Luckily, there is an alternative approach, as I report in my response to this message. Specifically, we can take advantage of the fact that the state for the root node of the tree computed during Felsenstein's independent contrasts algorithm also corresponds to the MLE for that node. We can thus go through all the nodes of the tree, each time re-rooting the tree at that node and computing its MLE ancestral character estimate by obtaining the PIC state for this new "root."
Here's how we might do this for phylogeny tree and data vector x:
> # first estimate alpha under the OU model
> fit<-fitContinuous(tree,x,model="OU")$Trait1
Fitting OU model:
> ou<-vector(); N<-length(tree$tip); M<-tree$Nnode
> # transform the tree by alpha
> outree<-ouTree(tree,fit$alpha)
> # compute the PIC "root" state for each internal node
> for(i in 1:M+N){
ou[i-N]<-ace(x,multi2di(root(outree,node=i)),
method="pic")$ace[1]
names(ou)[i-N]<-i
}
Simple as that. This will give the same estimates as ace(...,method="ML") (or, for that matter, phytools::anc.ML) for conditions in which the likelihood can be maximized by R.
Note that for some reason that is not completely clear to me root(...,resolve.root=TRUE) and multi2di(root(...)) don't yield the same tree & branch lengths - the latter being what we need in this case.
Hi Liam,
ReplyDeletethanks for this post! My comment is probably a bit late, but I wonder: how (and why!) does ouTree work at all?
I understand that it is transforming branch lengths so that the ML estimate for standard BM will return ancestors under the OU model - but how is that possible?
I can imaging that for a character that grows, say, quadratically or exponentially over time, you modify branches by lengthening them appropriately, and I can somehow imagine it for Pagels lambda etc., but I don't get it for OU.
I looked into the code and it does something like scale each branch exponentially according to its length and height in the tree and according to alpha. But I wonder: there is no theta involved in that formula.
I know that the tip values under OU are distributed according to a normal distribution, so in that sense it is similar to a BM model. But while the mean of that Normal under a BM is simply the value of the root, the mean of the Normal under an OU model is dependent on theta.
So I wonder - how can scaling the branches and using ML inference be enough to infer ancestral states under OU - since (I believe) this scaling can only influence the variation, but not the mean? Is there a manuscript somewhere that explains the math behind this?
I'm sure I'm missing something very basic here, but I can't figure it out.
It seems that transforming the tree and then using standard ML under BM is a standard way to find a ML estimate under more complex models - but how does this magic work?
Thanks a lot!
Best,
Christina