In a recent R-SIG-phylo thread, a R phylogenetics user inquired about how to obtain ancestral state estimates for internal nodes based on an Ornstein-Uhlenbeck model for trait evolution on the tree. They pointed out that, in principle, it should be possible to do this using ace(...,method="GLS",corStruct=corMartins(...)); unfortunately this was not producing sensible ancestral character estimates (and I wasn't able to get it to work at all, for some unknown reason).
Without judging the wisdom of this general endeavor (which I'll discuss a little at the end of this post), I provided a relatively straightforward way to do this using ace(...,method="ML") combined with geiger::fitContinuous and geiger::ouTree. Basically, we first fit our OU model using fitContinuous(...,model="OU"), then we transform the branch lengths of our tree using ouTree(...), finally we estimate ancestral states using ace in the 'ape' package or phytools::anc.ML.
Just to demo this, I'll simulate and estimate below:
> # load the required libraries
> library(phytools) # also loads ape
> # set the conditions for simulation
> N<-50; alpha<-1
> # simulate our tree
> # simulate data, including internal nodes
> # under a single-optimum OU model
> # fit using geiger::fitContinuous
Fitting OU model:
> # estimate ancestral states
51 52 53 54
0.012733299 0.012733173 0.012733128 ...
> # or ou<-anc.ML(ouTree(tree,fit$alpha),x[1:N])
> # check how we did
> plot(x[names(ou$ace)],ou$ace,xlab="true states",ylab="estimated states")
> plot(x[names(bm$ace)],bm$ace,xlab="true states",ylab="estimated states")
I'd like to estimate my ancestral states, which evolve under OU model.
But before using the tool you presented in this post, I need to estimate the "alpha" parameter.
But fitContinuous (fit<-fitContinuous(tr,dat,model="OU")) and anc.ML (fit<-anc.ML(tr,dat,model="OU")) give widely different values.
fitted ‘OU’ model parameters:
alpha = 0.467411
sigsq = 974.957102
z0 = 196.516937
log-likelihood = -366.988318
Fitted model parameters & likelihood:
sigma^2 alpha logLik
76.2096 0.022422 -642.7123
Which one should I use ? Why such a difference ?
Tks for this blog,