Hi Emmanuel, Many thanks for your suggestions. I think a group of us at Berkeley would be very interested in working out different ways to simulate these processes and then share our results. We'll keep you posted. I think starting with rcoal and compute.brlen would be a great jumping off place.
Thanks again, Andy On Mon, Apr 11, 2011 at 4:36 AM, Emmanuel Paradis <emmanuel.para...@ird.fr>wrote: > 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/ > [[alternative HTML version deleted]] _______________________________________________ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
