Hello Dave,
from my part, yes, this does exactly what I wanted to do. As far as I
understood, it is also what Rob was suggesting.
Thanks to everyone for helping with this.
Cheers,
Antigoni
On 10/13/2011 07:02 AM, David Bapst wrote:
Hello Rob and Antigoni,
Would the following meet your purposes?
Cheers,
-Dave Bapst, UChicago
library(ape)
tree<-rtree(30) #you can use your own tree here or rcoal() or whatever
Nchanges<-10 #number of shifts at nodes
nodeShifts<-sample((2:Nnode(tree))+Ntip(tree),Nchanges)
DescList<-prop.part(tree)[nodeShifts-Ntip(tree)]
tipChar<-rep(T,Ntip(tree))
for(i in 1:length(DescList)){
tipChar[DescList[[i]]]<-!tipChar[DescList[[i]]]
}
tipChar<-as.numeric(tipChar)
plot(tree,show.tip.label=F)
tiplabels(pch=21,bg=tipChar,cex=2)
nodeShifts #which nodes the changes occurred at
tipChar #characters at tips
On Wed, Oct 12, 2011 at 10:33 PM, Rob Lanfear<[email protected]> wrote:
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
--
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
