ancestral.pml gives NaN values

[kbhoehn]

Hi Klaus,

Thanks for maintaining this package. I've encountered issues in which ancestral.pml returns NaN values when there is data for every sequence at every site. This seems to be related to using di2mult before ancestral.pml, however I get similar issues with large numbers of sequences. Please see the attached code and file to reproduce the error. I'm using phangorn 2.12.1 and ape 5.8. Further, if I use return="phyDat" only 4 sites are returned.

Best,
-Ken
ancestral_pml_error.zip

[KlausVigo]

Hi @kbhoehn,

That seems to be caused by some numerical over / underflow. In pml I take care of it by scaling the likelihoods. In a few cases I use a simpler approach. Seems to look like I need to take more care and also introduce scaling with ancestral reconstruction, especially for larger phylogenies.
In your data as there are many duplicated sequences, in fact there are only 3 unique sequences and and we can build the reconstruction from only these sequences without losing any information. In such cases using a multifurcating tree avoids this problem and the likelihood etc. should be exactly the same.

Kind regards,
Klaus

PS:
It is not a good idea to change objects from a pml object directly, it is safer to use

fit <- update(fit, tree=tree)

as the tree is expected to be in "postorder".

PPS:
phangorn now contains also a joint reconstruction hidden in the namespace phangorn:::joint_pml and used inside anc_pml, whereas return="phyDat" assigns the state with the highest posterior. joint_pml` is limited for now, it does currently not allow site heterogeneity (gamma model, invariant sites etc.).