Hi Yan,
multi2di() calls rtree() if random = TRUE. As you rightly guessed, the
algorithm used by this latter function is (described in APER2, p. 313):
1. Draw randomly an integer a on the interval [1, n − 1]. Set b = n − a.
2. If a > 1, apply (recursively) step 1 after substituting n by a.
3. Repeat step 2 with b in place of a.
4. Assign randomly the n tip labels to the tips.
This is indeed results in unequal frequencies of the possible topologies. If
step 1 "n − 1" is replaced by "floor(n/2)", then all topologies are generated
with equal probability.
Practically, I see two possibilities for you. If pratical, you generate a list
with all possible topologies, for instance with phangorn's function allTrees:
R> library(phangorn)
R> TR <- allTrees(4, rooted = TRUE)
R> TR
15 phylogenetic trees
Then you just draw randomly some integers between 1 and 15 and use them as
indices:
R> s <- sample.int(length(TR), size = 1e5, replace = TRUE)
R> TR[s]
100000 phylogenetic trees
If you want, you can play with the option 'prob' of sample.int().
The other solution is to make copies of rtree() and multi2di() in your
workspace: only rtree() needs to be modified and only where sample.int() is
called in the recursive function foo. You also need to copy multi2di() even
unchanged because of the package's namespace.
HTH
Best,
Emmanuel
----- Le 27 Juin 20, à 4:46, Yan Wong y...@pixie.org.uk a écrit :
> I’m trying to figure out how to randomly resolve polytomies such that there is
> an equal probability of any topology being generated. I thought that the ape
> function “multi2di” did this, but when I have tried it repeatedly with a
> 4-tomy, multi2di seems to produce the 3 balanced trees [((a,b),(c,d)) ;
> ((a,c),(b,d)); ((a,d),(b,c))] twice as often as the 12 possible unbalanced
> dichotomous 4-tip rooted topologies. The R code I’ve used to produce the
> sample
> topologies is something like this:
>
> do.call(c, lapply(1:100000, function(x)
> multi2di(starTree(c('a','b','c','d')))))
>
> Firstly, is this expected, or am I doing something wrong (if expected, it
> would
> be useful to note this in the docs)? Secondly, is there an function somewhere
> that *will* break polytomies to produce equiprobable topologies? If not,
> thirdly, is there an algorithm that will do this? I think the standard
> “repeatedly pick 2 random edges from the polytomy node and pair them off”
> results in the non-equiprobable distribution that I find using multi2di. I
> think I’ve found a similar problem with the Dendropy algorithm, which does
> claim to result in equiprobable topologies, and have posted to their mailing
> list in case I’m misunderstanding something.
>
> Cheers
>
> Yan Wong
> Big Data Institute, Oxford University
> _______________________________________________
> R-sig-phylo mailing list - R-sig-phylo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
