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
>
--
David Bapst
Dept of Geophysical Sciences
University of Chicago
5734 S. Ellis
Chicago, IL 60637
http://home.uchicago.edu/~dwbapst/
_______________________________________________
R-sig-phylo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
