Andy,
Andy Rominger wrote on 09/04/2011 08:37:
Hello,
I'm working with a friend on a problem where we'd like to simulate
many possible trees that all have the same topology, but differ in branch
lengths (or equivalently in node ages). We have a known age for the root,
and all tips correspond to extant taxa. So the start and end times are
known.
The code to do this is in rcoal() except that a random topology is
generated. It seems several users are interested in assigning random
branching times to a fixed topology, so I'll add this. compute.brlen()
can generate random branch lengths on a fixed topology. What you want is
a function that does the opposite operation of branching.times(). I can
see three possible ways to implement this:
1) add a new option to compute.brlen()
2) add a new option to rcoal()
3) write a new function, eg, compute.brtime()
BTW, rcoal() can be used with fixed branching intervals (ie, computed
beforehand) where a random coalescent topology is generated; so option 2
makes sense to me. Any suggestion?
I'm wondering if a function to produce random branch lengths under
topological constraints already exists? If not, I'm trying to dredge up my
foggy "knowledge" of stochastic processes to write the simulation myself.
I'm thinking that maybe I remember under a Yule process sojourn times are
distributed exponential so event times should be uniformly distributed???
But then I got pointed toward "Yule priors" which I still
don't satisfactorily understand (but maybe eventually will) after a few web
searches.
Once you'll have the above tool (soon I hope), then you'll be able to
explore this issue more freely. Be careful to distinguish branching
times and branching intervals. (I'm just seeing that the help page of
rcoal is not clear on this point: I'll clarify this.) If you simulate
the latter, you get the former with cumsum(). ape has an internal
function that computes the theoretical CDF of branching times under any
birth-death process. eg:
CDF.birth.death(.1, 0, Tmax = 10, x = 0:100/10, case = 1)
will return the cumulative probability of the branching times 0, 0.1,
0.2, ..., 10 for lambda = 0.1 and mu = 0 (thus a Yule process). From
this you can simulate a set of random branching times with the inverse
method (generate runif(n) and back transform from the CDF). The
literature reference is in ?bd.time. Of course, you can use another
formula and compare the results with standard plots in R. As a reminder,
rcoal() generates coalescent intervals with:
2 * rexp(n - 1)/(as.double(n:2) * as.double((n - 1):1))
You may also have a look at the package TreeSim which has algorithms
contributed by Tanja Stadler.
Cheers,
Emmanuel
Any insights would be greatly appreciated. Thanks in advance--
Andy
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
Emmanuel Paradis
IRD, Jakarta, Indonesia
http://ape.mpl.ird.fr/
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
