Ancestral state reconstructions are incorrect when pis are not equal

[arklumpus]

Hi, I think there might be a bug in ancestral state estimates made using make.simmap, when the pis are not equal.

Consider the following example:

library("phytools")

tree <- read.tree(text="(A:1,B:1);");
data <- setNames(c(0, 1), c("A", "B"));
simmap <- make.simmap(tree = tree, x = data, model="ER", nsim=10000, pi=c(0.7, 0.3));
simmap_summary <- summary(simmap);

simmap_summary$ace

This is a very simple tree (just two branches with length 1) with a binary character. Since the branches have the same lengths and the model has equal rates, the likelihood at the root node is identical for the two states. As a result, the posterior reconstruction for the ancestral state of the root should be equal to the prior.

However, the output of the previous code is something like:

       0      1
3 0.8475 0.1525
A 1.0000 0.0000
B 0.0000 1.0000

This might be because lines 248-249 in make.simmap.R assign the root node state by multiplying the likelihood L by the prior pi; however, the likelihoods at the root node have already been multiplied by pi in line 105 of fitMk.R.

Indeed, , which is quite similar to the value I got using the example above (likely within the variability of the simulation process).

If the tree has more than two tips, I think this affects the estimate at all internal nodes, since those incorrect root posteriors are used to draw the initial states for the simulations.

Of course, there is no issue if the pis are equal, but I think the same problem applies even if the they are estimated or fitzjohn, instead of being manually fixed.

[liamrevell]

Thank you. I think you're right. Thanks so much for bringing this to my attention. I am doublechecking to make sure I understand that you are correct & then I will fix & push to GitHub & CRAN.

[liamrevell]

Fixed on GitHub now. Will also update on CRAN soon.

[arklumpus]

Great, thanks!