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/. 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
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 10/11/2011 6:18 PM, Antigoni Kaliontzopoulou wrote:
Hello everyone,
I am studying a binomial trait across a phylogeny and I'm trying to
simulate specific numbers of evolutionary shifts (1, 5, 10 shifts)
across random trees.
Does anyone know how this can be done?
Any help is greatly appreciated
Cheers,
Antigoni
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
