Re: [R-sig-phylo] Simulating binomial trait shifts on a phylogeny
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 rob.lanf...@gmail.com 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 antig...@mail.icav.up.pt
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: antig...@mail.icav.up.pt
ant...@gmail.com
Re: [R-sig-phylo] Simulating binomial trait shifts on a phylogeny
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 Lanfearrob.lanf...@gmail.com 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 Kaliontzopoulouantig...@mail.icav.up.pt
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
Re: [R-sig-phylo] Simulating binomial trait shifts on a phylogeny
Hello again everyone,
forwarded is some code provided by Dave Bapst for simulating a certain
number of trait shifts in nodes AND tracking those changes as you go.
This code will simulate Nchanges shifts and give both tip (tipChar)
and node (intNodeChar) states at the end of the run.
Dave correctly suggested this might be useful for other people, so I am
forwarding his code with my apologies for discussing this with him offline.
Dave, thanks again for the help. I think this can be really helpful for
some simulation purposes.
Cheers,
Antigoni
On 10/13/2011 02:53 PM, David Bapst wrote:
Hello Antigoni,
I think most people might actually want to know that, but I didn't
code it to make it easy to figure it out. So, I can understand why you
had trouble getting my simple brutish kludge to do it. The following
code should do what you are looking for. If you find this helpful,
please forward this email to R-Sig-Phylo for posterity.
-Dave
#a different version which will track the character changes
library(phangorn)
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-Descendants(tree,nodeShifts,type=all)
nodeChar-rep(T,Ntip(tree)+Nnode(tree))
for(i in 1:length(DescList)){
nodeChar[DescList[[i]]]-!nodeChar[DescList[[i]]]
nodeChar[nodeShifts[i]]-!nodeChar[nodeShifts[i]]
}
nodeChar-as.numeric(nodeChar)
tipChar-nodeChar[1:Ntip(tree)]
intNodeChar-nodeChar[-(1:Ntip(tree))]
plot(tree,show.tip.label=F)
tiplabels(pch=21,bg=tipChar,cex=2)
nodelabels(pch=21,bg=intNodeChar,cex=2)
nodeShifts #which nodes the changes occurred at
tipChar #characters at tips
nodeChar#character at each node, including tips
intNodeChar #character at each internal node
On Thu, Oct 13, 2011 at 1:50 PM, Antigoni Kaliontzopoulou
antig...@mail.icav.up.pt wrote:
Hi Dave,
another question (off-R-sig-phylo because I don't think it's particularly
interesting for people): do you by any chance also have any code that will
track the final node character states after simulating the changes?
I've been banging my head on that for some hours now, so I'd really
appreciate any help with that...
Thanks a lot!
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 Lanfearrob.lanf...@gmail.com
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
Kaliontzopoulouantig...@mail.icav.up.pt
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
Re: [R-sig-phylo] Simulating binomial trait shifts on a phylogeny
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/. 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: antig...@mail.icav.up.pt
ant...@gmail.com
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] Simulating binomial trait shifts on a phylogeny
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 antig...@mail.icav.up.pt
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: antig...@mail.icav.up.pt
ant...@gmail.com
__**_
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/**listinfo/r-sig-phylohttps://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
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] Simulating binomial trait shifts on a phylogeny
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
