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

Reply via email to