Addendum to estimating ancestral states under an OU model

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 σ² 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.