Hi Antigoni, I had been thinking about some of these kinds of simulations recently too. One challenge I encountered (which is similar to your problem here, if I've understood it) was how to make my simulations actually relevant.
Specifically, given a particular tree I needed to simulate data with N changes in a trait on the phylogeny, and I needed to make sure that what I simulated was something I would have recovered using the simulated data at the tips. This won't always be the case for Liam's method (if I've understood it correctly), because even if I fix the tree (trivial) I could end up with simulated data (e.g. multiple changes on a single branch) which I would not then reconstruct using the data from the tips. Note that this isn't a problem with Liam's method, but would be a problem if I wanted to use it for my purposes. I gave up in my attempt to solve this, but this thread, and Liam's approach, made me wonder if this problem could be solved using a related approach. In particular, one could use a fixed tree, and change the loop in Liam's method to randomise the tip labels, re-do the ancestral state reconstruction, and keep going until you get a reconstruction with the desired number of changes. That would make sure that you only simulate datasets which you would have recovered with extant data, and it also has the advantage (for me at least) of holding the topology and the number of extant species with each trait value constant across replicates. I'd be interested to know if anyone thinks this is a reasonable approach? Cheers, Rob On 13 October 2011 02:20, Antigoni Kaliontzopoulou <[email protected] > wrote: > Hi Liam, > > thank you very much for this. Yes, the function works beautifully for what > we wanted to do and the code you sent greatly facilitates implemantation. > I was, however, wondering whether you can somehow "limit" the shifts to > happen in nodes, to avoid character reversals across a single branch. > > Again, thanks a lot for the help. > > Cheers, > Antigoni > > > On 10/11/2011 08:34 PM, Liam J. Revell wrote: > >> Hi Antigoni. >> >> I have a function in the new version of "phytools" called sim.history() >> that might help with this. This version is not available yet on CRAN, but >> it can be downloaded from my webpage: http://anolis.oeb.harvard.edu/** >> ~liam/R-phylogenetics/<http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/>. >> The function generates stochastic histories according to some rate matrix, >> and stores the histories along the states at all terminal nodes. >> >> This code might be helpful in your task. It simulates random trees, for a >> given number of changes it computes the instantaneous rate matrix that will >> produce that number of changes, and then it simulates (repeatedly) under >> that rate until the desired number of changes is obtained. The code >> fragment uses functions in "ape", "geiger", and "phytools": >> >> nchanges<-5 # desired number of changes >> ntrees<-100 # desired number of trees >> ntaxa<-100 # desired number of taxa >> require(phytools); require(geiger) # required packages >> mtrees<-list() >> for(i in 1:ntrees){ >> # simulate a stochastic B-D tree using birthdeath.tree >> >> tree<-drop.tip(birthdeath.**tree(b=1,d=0,taxa.stop=(ntaxa+**1)),as.character(ntaxa+1)) >> >> # compute the rate that results in (on average) nchanges changes >> q<-nchanges/sum(tree$edge.**length) >> # simulate a stochastic history >> mtree<-sim.history(tree,Q=**matrix(c(-q,q,q,-q),2,2)) >> # loop until desired number of changes >> while((length(unlist(mtree$**maps))-nrow(mtree$edge))!=**nchanges) >> mtree<-sim.history(tree,Q=**matrix(c(-q,q,q,-q),2,2)) >> mtrees[[i]]<-mtree >> # for fun, plot the history >> plotSimmap(mtrees[[i]]) >> } >> class(mtrees)<-"multiPhylo" >> # tip states for the kth tree are in mtrees[[k]]$states >> >> Please let us know if this is useful or if you have any questions or >> difficulties implementing. >> >> All the best, Liam >> >> > -- > Antigoni Kaliontzopoulou > > CIBIO, Centro de Investigacao em Biodiversidade e Recursos Geneticos > Campus Agrario de Vairao, 4485-661 Vairao > PORTUGAL > Department of Ecology, Evolution, and Organismal Biology > Iowa State University, Ames, > Iowa 50011, USA > > tel: +351 91 3086188 > mail to: [email protected] > [email protected] > > ______________________________**_________________ > R-sig-phylo mailing list > [email protected] > https://stat.ethz.ch/mailman/**listinfo/r-sig-phylo<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> > -- Rob Lanfear Postdoc, Centre for Macroevolution and Macroecology, Research School of Biology, Australian National University Tel: +61 2 6125 7270 www.robertlanfear.com [[alternative HTML version deleted]] _______________________________________________ R-sig-phylo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
