Re: [R-sig-phylo] selecting a set of incongruent trees from a posterior distribution
Hi Jesse, As Eduardo says, if in fact you want to see how different trees are from a "consensus", something that you could try is find the maximum-clade-credibility (MCC) tree (you can do this with treeAnnotator from BEAST). This will be a fully bifurcating tree and is essentially the tree in the posterior distribution that maximizes the overall clade posterior probabilities (e.g. sum(posterior)). Once you have the MCC tree, you can use the Robinson-Foulds distance metric (there is a function in phytools) to compare all the topologies in your posterior distribution. The trees with the lowest values will be the most incongruent with respect to the MCC tree. Just remember to exclude the burnin. Cheers, Santiago On Wed, Jul 26, 2017 at 6:01 AM wrote: > Send R-sig-phylo mailing list submissions to > r-sig-phylo@r-project.org > > To subscribe or unsubscribe via the World Wide Web, visit > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > or, via email, send a message with subject or body 'help' to > r-sig-phylo-requ...@r-project.org > > You can reach the person managing the list at > r-sig-phylo-ow...@r-project.org > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of R-sig-phylo digest..." > > > Today's Topics: > >1. selecting a set of incongruent trees from a posterior > distribution (Jesse Delia) >2. Re: selecting a set of incongruent trees from a posterior > distribution (Eduardo Ascarrunz) > > > -- > > Message: 1 > Date: Tue, 25 Jul 2017 16:16:30 -0400 > From: Jesse Delia > To: R-sig-phylo@r-project.org > Subject: [R-sig-phylo] selecting a set of incongruent trees from a > posterior distribution > Message-ID: > < > ca+lom6ehmv6c2v_lso8b2sdqqvoqwfcrl56-4jhcmicvjyv...@mail.gmail.com> > Content-Type: text/plain; charset="UTF-8" > > Dear list, > > Does anyone know of a function or script that could select a set of the > most incongruent trees from a list of trees? Maybe I missed a post, but > haven't found anything. > > I running mixed models, which take a lot of computational space on my lap > top. In effort to account for phylogenetic uncertainty, without having to > run 100s of chains, I figured maybe i could run models across a (much) > shorter list that accounts for more diversity in tree shape observed within > the posterior distribution. Not sure if this makes sense, and/or is > extremely complicated? > > Thanks for you time! > > Jesse > > [[alternative HTML version deleted]] > > > > -- > > Message: 2 > Date: Wed, 26 Jul 2017 11:33:12 +0200 > From: Eduardo Ascarrunz > To: Jesse Delia > Cc: R Sig Phylo Listserv > Subject: Re: [R-sig-phylo] selecting a set of incongruent trees from a > posterior distribution > Message-ID: > < > cakxhxpxj_8ucosjt-bzpzvsn+5quet9vxnd4ssd4bqaq-fp...@mail.gmail.com> > Content-Type: text/plain; charset="UTF-8" > > Hi Jesse, > > Do you want to select the trees that are most incongruent with some sort of > average of the tree distribution? The `treespace`and `phytools` packages > have functions that compute average trees with different metrics, and they > also allow you to measure the distance between each tree and the average > tree. > > Cheers, > > Eduardo > > 2017-07-25 22:16 GMT+02:00 Jesse Delia : > > > Dear list, > > > > Does anyone know of a function or script that could select a set of the > > most incongruent trees from a list of trees? Maybe I missed a post, but > > haven't found anything. > > > > I running mixed models, which take a lot of computational space on my lap > > top. In effort to account for phylogenetic uncertainty, without having to > > run 100s of chains, I figured maybe i could run models across a (much) > > shorter list that accounts for more diversity in tree shape observed > within > > the posterior distribution. Not sure if this makes sense, and/or is > > extremely complicated? > > > > Thanks for you time! > > > > Jesse > > > > [[alternative HTML version deleted]] > > > > ___ > > R-sig-phylo mailing list - R-sig-phylo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > > Searchable archive at http://www.mail-archive.com/r- > > sig-ph...@r-project.org/ > > > > [[alternative HTML version deleted]] > > > > -- > > Subject: Digest Footer > > ___ > R-sig-phylo mailing list > R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > > -- > > End of R-sig-phylo Digest, Vol 114, Issue 10 > > -- == Santiago Sanchez-Ramirez, PhD Postdoctoral Associate Ecology and Evolutionary Biology University of Toronto == [[alternativ
Re: [R-sig-phylo] need finite 'xlim' values error with plot.phylo
Thanks, Liam! On Tue, Jun 13, 2017 at 12:48 PM Liam Revell wrote: > This is a bug that is related to how the function optimizes xlim to permit > space for the tip labels. By not showing tip labels, there is no > optimization, and thus no error. > > > > This doesn’t seem to create an issue for plotTree in phytools, which has > many of the same options (and some additional ones) as plot.phylo. > > > > -- > Liam J. Revell, Associate Professor of Biology > University of Massachusetts Boston > web: http://faculty.umb.edu/liam.revell > email: liam.rev...@umb.edu > > Sent from my Windows 10 phone > > > > *From: *Santiago Sánchez > *Sent: *Tuesday, June 13, 2017 11:35 AM > *To: *r-sig-phylo@r-project.org > *Subject: *[R-sig-phylo] need finite 'xlim' values error with plot.phylo > > > Hello, > > I was wondering if there could be a bug in the plot.phylo code. I'm trying > to plot a tree which is scaled in yr units (less than 1 million years). > However I'm getting this error: > > Error in plot.window(...) : need finite 'xlim' values > > I tried rescaling to kyr and it plots well. Here is a sample test code: > > library(ape) > rtr <- rtree(50) > rtr$edge.length <- rtr$edge.length*10 > plot(rtr) > rtr2 <- rtr > rtr2$edge.length <- rtr2$edge.length/1000 > plot(rtr2) > > Something interesting is that it will plot if I omit tip labels: > > plot(rtr, show.tip.label=F) > > Cheers, > Santiago > -- > == > Santiago Sanchez-Ramirez, PhD > Postdoctoral Associate > Ecology and Evolutionary Biology > University of Toronto > == > > [[alternative HTML version deleted]] > > ___ > R-sig-phylo mailing list - R-sig-phylo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > Searchable archive at > http://www.mail-archive.com/r-sig-phylo@r-project.org/ > -- == Santiago Sanchez-Ramirez, PhD Postdoctoral Associate Ecology and Evolutionary Biology University of Toronto == [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] need finite 'xlim' values error with plot.phylo
Hello, I was wondering if there could be a bug in the plot.phylo code. I'm trying to plot a tree which is scaled in yr units (less than 1 million years). However I'm getting this error: Error in plot.window(...) : need finite 'xlim' values I tried rescaling to kyr and it plots well. Here is a sample test code: library(ape) rtr <- rtree(50) rtr$edge.length <- rtr$edge.length*10 plot(rtr) rtr2 <- rtr rtr2$edge.length <- rtr2$edge.length/1000 plot(rtr2) Something interesting is that it will plot if I omit tip labels: plot(rtr, show.tip.label=F) Cheers, Santiago -- == Santiago Sanchez-Ramirez, PhD Postdoctoral Associate Ecology and Evolutionary Biology University of Toronto == [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] getBtimes vs. branching.times
Dear Lev,
I simple way to do what you request would be to follow this example with
your data:
> library(ape)
# here I'm creating a random tree with branch lengths
> treeBL <- rtree(10)
> treeBL$edge.length
[1] 0.59962263 0.76458987 0.48171189 0.26387301 0.76967429 0.42754456
[7] 0.30173612 0.10535179 0.01769042 0.49671428 0.61575421 0.91297819
[13] 0.52992782 0.68690369 0.35537999 0.56974770 0.97019530 0.49849792
# this would be the tree with the same topology but without branch lengths
> treeNN <- treeBL
> treeNN$edge.length <- NULL
> treeNN
Phylogenetic tree with 10 tips and 9 internal nodes.
Tip labels:
t1, t9, t4, t3, t7, t10, ...
Rooted; no branch lengths.
# create node labels
> treeNN$Nnode
[1] 9
> treeNN$node.label <- paste("n",1:treeNN$Nnode,sep='')
> treeNN$node.label
[1] "n1" "n2" "n3" "n4" "n5" "n6" "n7" "n8" "n9"
> treeNN
Phylogenetic tree with 10 tips and 9 internal nodes.
Tip labels:
t1, t9, t4, t3, t7, t10, ...
Node labels:
n1, n2, n3, n4, n5, n6, ...
Rooted; no branch lengths.
# check if the internal structure of both trees is identical
> identical(treeNN$edge,treeBL$edge)
[1] TRUE
> if (identical(treeNN$edge,treeBL$edge)){
+ newtree <- treeBL
+ newtree$node.label <- treeNN$node.label
}
> newtree
Phylogenetic tree with 10 tips and 9 internal nodes.
Tip labels:
t1, t9, t4, t3, t7, t10, ...
Node labels:
n1, n2, n3, n4, n5, n6, ...
Rooted; includes branch lengths.
Hope this helps.
Cheers,
Santiago
--
Santiago Sánchez-Ramírez
Environmental Genomics Group
Max Planck Institute for Evolutionary Biology
August-Thienemann-Str. 2
24306 Plön
Germany
Website: https://sites.google.com/site/santiagosnchezrmirez/
>
> Message: 1
> Date: Wed, 13 Jan 2016 23:39:21 +
> From: "Yampolsky, Lev"
> To: Emmanuel Paradis , "Liam J. Revell"
> , Dan Rabosky
> Cc: "r-sig-phylo@r-project.org"
> Subject: Re: [R-sig-phylo] getBtimes vs. branching.times
> Message-ID:
> Content-Type: text/plain; charset="utf-8"
>
> Dear Colleagues,
>
> Thank you very much for your help, it?s all clear now. (I got way more
> than I asked for, for example, my trees include only extant taxa). And
> yes, branching.times() and getBtimes() return exact same numbers of
> course, I just didn?t understand ow indexing works.
>
> Next question, a very simple one again.
>
> I have two newick trees, identical topology. In one I have branch lengths,
> in the other I have node names. I would like a tree with both. How do I do
> this?
>
> Thanks in advance!
>
> --
> Lev Yampolsky
>
> Professor
> Department of Biological Sciences
> East Tennessee State University
> Box 70703
> Johnson City TN 37614-1710
> Cell 423-676-7489
> Office/lab 423-439-4359
> Fax423-439-5958
>
>
>
>
>
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] simulate Ne changes on a given phylogeny
Hi Jacob, There is a package named OutBreakTools ( https://cran.r-project.org/web/packages/OutbreakTools/index.html) that includes a function named "read.annotated.nexus". It will retrieve all metadata from trees produced by BEAST/TreeAnnotator. I found particularly useful that you can also import "raw" posterior trees with their metadata. I've tried it on species trees from *BEAST and it works well. In contrast to "read.beast" from phyloch, "read.annotated.nexus" will associate metadata to edges (e.g. edge number), rather than nodes. Cheers, Santiago -- Santiago Sánchez-Ramírez, PhD Department of Ecology and Evolutionary Biology, University of Toronto Department of Natural History (Mycology), Royal Ontario Museum 100 Queen's Park Toronto, ON M5S 2C6 Canada Other email: santiago.sanc...@mail.utoronto.ca Tel. 416-586-8025 Website: https://sites.google.com/site/santiagosnchezrmirez/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] simulate Ne changes on a given phylogeny
Hi Jacob, You might want to check Joseph Heled's python package biopy (https://code.google.com/p/biopy), I know it has some simulation capabilities. Also check Z. Yang's program MCMCcoal, I know some people have used it for simulation analyses in BPP. Cheers, Santiago Sent from my iPhone ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Add a random outgroup lineage to a phylo object
Hello everyone,
I thought someone might find this useful. I went on to write a modest
function that adds a random outgroup lineage to a phylo object. It accepts
both ultrametric and standard trees with or without branch lenghts. In the
case of ultrametric trees, it will add an "ingroup" edge (leading to the
root) of length: 1/5*rootheigth, where rootheight is the oldest branching
time of the original tree. And in the case of standard trees with branch
lengths, the "outgroup" branch length will be as long as the longest branch
in the tree, and the ingroup edge will be 1/3 as long. You can also reroot
your tree internally by giving the node number of the new root (root.node).
Its possible that Liam Revell already has something similar, but I'm not
completely sure. In case you are feeling curious, I prepared a quick and
simple demo:
# The code:
add.outgroup <- function(phy, root.node=NULL){
require(ape)
if (is.null(phy$edge.length)){
tt <- rtree(2)
tt$edge.length <- NULL
tt$tip.label <- c("outgroup","drop")
phy$root.edge <- 1
if (!is.null(root.node)){
rphy <- root(phy, node=root.node)
ot <- bind.tree(rphy, tt, position=1)
} else {
ot <- bind.tree(phy, tt, position=1)
}
ot <- bind.tree(phy, tt, position=1)
res <- drop.tip(ot, "drop")
return(res)
} else if (is.ultrametric(phy)){
th <- as.numeric(sort(branching.times(phy), decreasing=T))[1]
re <- th/5
phy$root.edge <- re
tt <- rtree(2)
tt$edge.length <- c(0,0)
tt$tip.label <- c("outgroup","drop")
tt$root.edge <- th + re
ot <- bind.tree(phy, tt, position=re)
res <- drop.tip(ot, "drop")
return(res)
} else {
tl <- max(phy$edge.length)
re <- tl/3
tt <- rtree(2)
tt$edge.length <- c(0,0)
tt$tip.label <- c("outgroup","drop")
tt$root.edge <- tl
if (!is.null(root.node)){
rphy <- root(phy, node=root.node)
rphy$root.edge <- re
ot <- bind.tree(rphy, tt, position=re)
} else {
phy$root.edge <- re
ot <- bind.tree(phy, tt, position=re)
}
res <- drop.tip(ot, "drop")
return(res)
}
}
# Demo:
> require(ape)
> # for a regular tree with branch lengths
> phy <- rtree(20)
> nodes <- (length(phy$tip.label)+1):(length(phy$tip.label)+phy$Nnode)
> root.node <- sample(nodes, 1)
> rphy.w.out <- add.outgroup(phy=phy, root.node=root.node)
> is.rooted(rphy.w.out)
> # try it without rooting
> phy.w.out <- add.outgroup(phy=phy)
> is.rooted(phy.w.out)
> # tree with no branch lengths
> phynbl <- phy
> phynbl$edge.length <- NULL
> phynbl.w.out <- add.outgroup(phy=phynbl)
> is.rooted(phynbl.w.out)
> par(mfrow=c(1,4))
> plot(ladderize(phy.w.out,FALSE), main="original")
> add.scale.bar(y=19, x=max(phy.w.out$edge.length)/2)
> plot(ladderize(rphy.w.out,FALSE), main="rerooted")
> add.scale.bar(y=19, x=max(rphy.w.out$edge.length)/2)
> plot(ladderize(phynbl.w.out,FALSE), main="no branch lengths")
> # for an ultrametric tree
> phy <- pbtree(n=20)
> nodes <- (length(phy$tip.label)+1):(length(phy$tip.label)+phy$Nnode)
> root.node <- sample(nodes, 1)
> phy.w.out <- add.outgroup(phy=phy, root.node=root.node)
> is.rooted(phy.w.out)
> plot(ladderize(phy.w.out,FALSE), main="ultrametric")
> axisPhylo()
Cheers,
Santiago
--
Santiago Sánchez-Ramírez
Department of Ecology and Evolutionary Biology, University of Toronto
Department of Natural History (Mycology), Royal Ontario Museum
100 Queen's Park
Toronto, ON
M5S 2C6
Canada
Other email: santiago.sanc...@mail.utoronto.ca
Tel. 416-586-8025
Website: https://sites.google.com/site/santiagosnchezrmirez/
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Help interpreting diversitree:MuSSe results
Hi Tim, I've used BiSSE in a similar context as you (e.g. in a tree with 62 taxa), and in my case allowing extinction rates to vary had a significantly higher AIC compared to simpler models such as extinction = 0 (Sanchez-Ramirez et al. 2015 Journal of Biogeography). In BiSSE type analyses (including MuSSE) rate estimation depends on the tree topology, the tree size, the branch lengths, and trait distribution along the tips. If if extinction = 0 is always favor, it probably means that in other models where extinction is allowed to vary, the estimate is not significantly different from 0. The other situation is that the power to detect any signature in the data also depends on the number of taxa (see FitzJohn 2012). Hope this helps, Santiago Message: 1 > Date: Sat, 24 Jan 2015 21:04:18 + > From: Tim Astrop > To: r-sig-phylo@r-project.org > Subject: [R-sig-phylo] Help interpreting diversitree:MuSSe results > Message-ID: > < > canw8cfbdtdz30kuf7uy16jposhdcq6y4jrv4yjyyc7wp86t...@mail.gmail.com> > Content-Type: text/plain; charset="UTF-8" > > Hello all, > > I'm playing with the MuSSe part of the Diversitree package looking at how a > 3 state character might affect diversity in a tree containing 56 taxa. > > I have some biological info to inform constraints on character transitions > and am trying to be holistic be looking at a many models as I think are > viable. > > My question is, why does a model where all extinction rates = 0 always > outperform (in terms of AIC) others, even across a sub-sample of posterior > trees? I know estimating extinction rates is considered less reliable than > estimating speciation rates. Is this a problem specific to me or one in > general? > > I hope this makes sense! > > Cheers > > *Tim Astrop PhD <http://timastrop.wordpress.com/>* > Ammonoid Palaeobiology Lab <http://aplbath.wordpress.com/> > Department of Biology & Biochemistry > University of Bath > > *** > -- Santiago Sánchez-Ramírez Department of Ecology and Evolutionary Biology, University of Toronto Department of Natural History (Mycology), Royal Ontario Museum 100 Queen's Park Toronto, ON M5S 2C6 Canada Other email: santiago.sanc...@mail.utoronto.ca Tel. 416-586-8025 Website: https://sites.google.com/site/santiagosnchezrmirez/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] profiles.plot legends
o solve.
>
> Cheers,
>
> Kev
>
>
> From: Nicholas Crouch [ncro...@uic.edu]
> Sent: 14 March 2014 16:46
> To: Arbuckle, Kevin
> Subject: Re: [R-sig-phylo] profiles.plot legends
>
> This thread might be useful as well
>
>
> http://www.talkstats.com/showthread.php/16691-greek-letter-and-equal-sign-in-legend-of-a-plot
>
> Nick
>
>
> On Fri, Mar 14, 2014 at 11:41 AM, Arbuckle, Kevin <
> k.arbuc...@liverpool.ac.uk<mailto:k.arbuc...@liverpool.ac.uk>> wrote:
> Thanks Nick, I am beginning to suspect that adding the legend separately
> is the way to go.
>
> Kev
>
>
> From: Nicholas Crouch [ncro...@uic.edu<mailto:ncro...@uic.edu>]
> Sent: 14 March 2014 14:59
> To: Arbuckle, Kevin
> Subject: Re: [R-sig-phylo] profiles.plot legends
>
> Kev,
>
> I'm pretty sure that would require going into the BiSSE code.
>
> Previously, I have just manually added the legend myself, then you can
> specify greek letters
>
> Nick
>
>
> On Fri, Mar 14, 2014 at 9:49 AM, Arbuckle, Kevin <
> k.arbuc...@liverpool.ac.uk<mailto:k.arbuc...@liverpool.ac.uk>> wrote:
> Hi all,
>
> I have a very simple question (I'm sure) that I can't quite seem to figure
> out. I have been trying to plot posterior distributions of parameter
> estimates from a BiSSE analysis using MCMC in diversitree. I have used the
> following code to compare the speciation rates between states (cdpostburn
> is the mcmc object), which generally does what I want.
>
> profiles.plot(cdpostburn[c("lambda0","lambda1")],col.line=c("blue",
> "red"),las=1,xlab="Speciation rate", legend="topright")
>
> The one thing is that the labels in the legend simply say 'lambda0' and
> 'lambda1', whereas I want them to be a little more informative. Again, I'm
> sure this is a very simple question but how do I specify the legend labels
> in the profiles.plot function?
>
> Thanks,
>
> Kev
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org R-sig-phylo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at
> http://www.mail-archive.com/r-sig-phylo@r-project.org/
>
>
>
> --
> Nicholas Crouch
> Graduate Student
> Department of Biological Sciences
> University of Illinois at Chicago
>
>
>
> --
> Nicholas Crouch
> Graduate Student
> Department of Biological Sciences
> University of Illinois at Chicago
>
> [[alternative HTML version deleted]]
>
>
>
> --
>
> Message: 6
> Date: Fri, 14 Mar 2014 15:05:30 -0400
> From: Jo Wolfe
> To: r-sig-phylo@r-project.org
> Subject: [R-sig-phylo] Outgroup with sim.char
> Message-ID:
> <
> cao6r+okf7xjspiey33obg4kgxno4spq+s_11-hdauz5ueu4...@mail.gmail.com>
> Content-Type: text/plain
>
> Hi R-ers,
>
> I am using the discrete argument in sim.char on a preset reference tree,
> and will eventually be using those simulated characters to re-reconstruct
> the tree in TNT, which requires a preset outgroup. In this case the
> outgroup must have all states as plesiomorphic (= state 1 in sim.char).
>
> Is the root initially set to all state 1, if root = TRUE? Otherwise, I
> assume I need to add an outgroup to the reference tree initially and use
> those states.
>
> Thanks!
> Jo
>
> --
> Jo Wolfe, Ph.D.
> Gerstner Scholar/Lerner-Gray Fellow
> Division of Invertebrate Zoology and Richard Gilder Graduate School
> America
--
Santiago Sánchez-Ramírez
Department of Ecology and Evolutionary Biology, University of Toronto
Department of Natural History (Mycology), Royal Ontario Museum
100 Queen's Park
Toronto, ON
M5S 2C6
Canada
Other email: santiago.sanc...@mail.utoronto.ca
Tel. 416-586-8025
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Fit BM, OU (1 theta), OU (2 theta) with the same function
Hi Christofer, Jeremy Beaulieu fixed this issue in OUwie v1.40. You might want to give it a try, although I'm not sure if it is already available on CRAN. Cheers, Santiago 2014-02-17 10:04 GMT-05:00 Santiago Sanchez : > Hi Christofer, > > Try with the package mvSLOUCH (available on CRAN). I ran my data with it > without a problem. The computation time might be a burden if you have many > large trees (~5 min per ~50 terminal tree and two optima). Also, if you > want to test the likelihood of having multiple regimes, it is easy to > convert for instance the phy$node.label vector that you would use in OUwie, > into the "regimes" vector used in mvSLOUCH. > > Jeremy Beaulieu also suggested log-transforming the trait values (in case > that you have many negative values) before running OUwie. > > I hope this helps. > > Cheers, > Santiago > > Santiago Sanchez-Ramirez > Ecology and Evolutionary Biology, University of Toronto > Natural History (Mycology), Royal Ontario Museum > 100 Queen's Park > Toronto, ON > M5S 2C6 > Canada > > > On Feb 17, 2014, at 6:00 AM, r-sig-phylo-requ...@r-project.org wrote: > > > > Send R-sig-phylo mailing list submissions to > >r-sig-phylo@r-project.org > > > > To subscribe or unsubscribe via the World Wide Web, visit > >https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > > or, via email, send a message with subject or body 'help' to > >r-sig-phylo-requ...@r-project.org > > > > You can reach the person managing the list at > >r-sig-phylo-ow...@r-project.org > > > > When replying, please edit your Subject line so it is more specific > > than "Re: Contents of R-sig-phylo digest..." > > > > > > Today's Topics: > > > > 1. Re: Fit BM, OU (1 theta),OU (2 theta) with the same function > > (Christofer Clemente) > > 2. Re: Fit BM, OU (1 theta),OU (2 theta) with the same function > > (Christofer Clemente) > > > > > > -- > > > > Message: 1 > > Date: Sun, 16 Feb 2014 11:00:38 + (UTC) > > From: Christofer Clemente > > To: r-sig-phylo@r-project.org > > Subject: Re: [R-sig-phylo] Fit BM, OU (1 theta),OU (2 theta) with the > >same function > > Message-ID: > > Content-Type: text/plain; charset=us-ascii > > > > > > > > > > > > Hi Brian and Santiago, > > > > I have been trying to test models of evolution using OUwie and i ran > into an > > identical problem to that Santiago describes. > > > > I am using a phylogeny (tree) of 111 species in ape's "phylo" format. > > my data frame (df) is of three columns: species, discrete trait, > continuous > > trait. > > > > bm1Output = OUwie(tree, df, model = "BM1") > > is successful > > > > ou1Output = OUwie(pruned, df1, model = "OU1") > > returns > > Error in svd(m) : infinite or missing values in 'x' > > > > ou1Output = OUwie(pruned, df1, model = "OU1", diagn=FALSE) > > similarly returns the error > > Error in svd(m) : infinite or missing values in 'x' > > > > Any suggestions greatly appreciated > > > > best > > > > Christofer > > > > > > > > -- > > > > Message: 2 > > Date: Mon, 17 Feb 2014 08:55:00 + (UTC) > > From: Christofer Clemente > > To: r-sig-phylo@r-project.org > > Subject: Re: [R-sig-phylo] Fit BM, OU (1 theta),OU (2 theta) with the > >same function > > Message-ID: > > Content-Type: text/plain; charset=us-ascii > > > > i should also note i am using the latest version of OUwie '1.39' > > > > > > > > -- > > > > ___ > > R-sig-phylo mailing list > > R-sig-phylo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > > > > > > End of R-sig-phylo Digest, Vol 73, Issue 11 > > *** > -- Santiago Sánchez-Ramírez Department of Ecology and Evolutionary Biology, University of Toronto Department of Natural History (Mycology), Royal Ontario Museum 100 Queen's Park Toronto, ON M5S 2C6 Canada Other email: santiago.sanc...@mail.utoronto.ca Tel. 416-586-8025 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Fit BM, OU (1 theta), OU (2 theta) with the same function
Hi Brian, Thank you for pointing these out. I've been playing around with OUwie and I ran into an error which I couldn't solve. This comes out when fitting any OU model implemented in OUwie, and it doesn't happen with BM models (they run fine): test <- OUwie(phy, data, model="OU1") Initializing... Finished. Begin thorough search... Finished. Summarizing results. Error in svd(m) : infinite or missing values in 'x' I would appreciate if you could give me some suggestions. Cheers, Santiago 2014-02-11 15:09 GMT-05:00 Brian O'Meara : > There are at least three packages to do this: OUCH, SLOUCH, and OUwie. > OUCH and OUwie are on CRAN; for SLOUCH, go to > http://freshpond.org/software/SLOUCH/ . > > Best, > Brian > > ___ > Brian O'Meara > Assistant Professor > Dept. of Ecology & Evolutionary Biology > U. of Tennessee, Knoxville > http://www.brianomeara.info > > Students wanted: Applications due Jan. 1, annually > Postdoc collaborators wanted: Check NIMBioS' website > Calendar: http://www.brianomeara.info/calendars/omeara > > > On Tue, Feb 11, 2014 at 3:05 PM, Santiago Sánchez < > santiago.snc...@gmail.com> wrote: > >> Hello list, >> >> Does someone know if there is a function out there that can fit (similar >> to >> fitContinuous() in geiger and compar.ou() in ape) OU models with multiple >> optima (theta) as well as a BM model. >> >> I'm not sure if fitContinuous() allows OU with multiple theta, and, as far >> as I know there is no "compar.bm" in ape. >> >> Any help is greatly appreciated. >> >> Cheers, >> Santiago >> >> -- >> Santiago Sánchez-Ramírez >> Department of Ecology and Evolutionary Biology, University of Toronto >> Department of Natural History (Mycology), Royal Ontario Museum >> 100 Queen's Park >> Toronto, ON >> M5S 2C6 >> Canada >> Other email: santiago.sanc...@mail.utoronto.ca >> Tel. 416-586-8025 >> >> [[alternative HTML version deleted]] >> >> >> ___ >> R-sig-phylo mailing list - R-sig-phylo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo >> Searchable archive at >> http://www.mail-archive.com/r-sig-phylo@r-project.org/ >> > > -- Santiago Sánchez-Ramírez Department of Ecology and Evolutionary Biology, University of Toronto Department of Natural History (Mycology), Royal Ontario Museum 100 Queen's Park Toronto, ON M5S 2C6 Canada Other email: santiago.sanc...@mail.utoronto.ca Tel. 416-586-8025 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Fit BM, OU (1 theta), OU (2 theta) with the same function
Hello list, Does someone know if there is a function out there that can fit (similar to fitContinuous() in geiger and compar.ou() in ape) OU models with multiple optima (theta) as well as a BM model. I'm not sure if fitContinuous() allows OU with multiple theta, and, as far as I know there is no "compar.bm" in ape. Any help is greatly appreciated. Cheers, Santiago -- Santiago Sánchez-Ramírez Department of Ecology and Evolutionary Biology, University of Toronto Department of Natural History (Mycology), Royal Ontario Museum 100 Queen's Park Toronto, ON M5S 2C6 Canada Other email: santiago.sanc...@mail.utoronto.ca Tel. 416-586-8025 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Error in read.beast()
To me it sounds like one of the BEAST annotation vectors is empty. Did you
inspect the file in FigTree and check if all annotations show up?
Cheers,
Santiago
On Tuesday, January 21, 2014, Mariana Vasconcellos
wrote:
> Thank you Santiago, but this didn't solve my problem. I have already tried
> this and it didn't work. I don't understand what the error: "Error in (i2 +
> 1):end : argument of length 0" means.
>
> --
> Mariana Mira Vasconcellos
> marian...@utexas.edu
> PhD candidate
> Ecology, Evolution & Behavior
> Integrative Biology
> The University of Texas at Austin
>
>
>
>
> On Jan 21, 2014, at 7:22 AM, Santiago Sánchez
> 'santiago.snc...@gmail.com');>>
> wrote:
>
> Hello Mariana,
>
> When I use read.beast() I usually specify a digit number.
> If most estimates are in millions of years I usually give digits=2, but if
> time frames are shorter, like thousands of years, you will have smaller
> fractions, in this case I give digits=10. Again, I'm not sure if this will
> solve the issue.
>
> > read.beast("tree", digits=2)
>
> Or
>
> > read.beast("tree", digits=10)
>
> Cheers,
> Santiago
>
>
> --
> Santiago Sánchez-Ramírez
> Department of Ecology and Evolutionary Biology, University of Toronto
> Department of Natural History (Mycology), Royal Ontario Museum
> 100 Queen's Park
> Toronto, ON
> M5S 2C6
> Canada
> Other email: santiago.sanc...@mail.utoronto.ca 'santiago.sanc...@mail.utoronto.ca');>
> Tel. 416-586-8025
>
>
>
--
Santiago Sánchez-Ramírez
Department of Ecology and Evolutionary Biology, University of Toronto
Department of Natural History (Mycology), Royal Ontario Museum
100 Queen's Park
Toronto, ON
M5S 2C6
Canada
Other email: santiago.sanc...@mail.utoronto.ca
Tel. 416-586-8025
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] Error in read.beast()
Hello Mariana,
When I use read.beast() I usually specify a digit number.
If most estimates are in millions of years I usually give digits=2, but if
time frames are shorter, like thousands of years, you will have smaller
fractions, in this case I give digits=10. Again, I'm not sure if this will
solve the issue.
> read.beast("tree", digits=2)
Or
> read.beast("tree", digits=10)
Cheers,
Santiago
--
Santiago Sánchez-Ramírez
Department of Ecology and Evolutionary Biology, University of Toronto
Department of Natural History (Mycology), Royal Ontario Museum
100 Queen's Park
Toronto, ON
M5S 2C6
Canada
Other email: santiago.sanc...@mail.utoronto.ca
Tel. 416-586-8025
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] Fwd: which.node function
-- Forwarded message --
From: *Santiago Sánchez*
Date: Monday, November 25, 2013
Subject: which.node function
To: "Liam J. Revell"
Hi Liam,
Great! What a waste of time! Anyway thank you pointing this out. I also
notice the error, and just in case, I'm posting a corrected version of the
code.
which.node <- function(phy, tips)
{
if (is.character(tips) == TRUE){
tips_n = vector("numeric",length(tips))
for (i in 1:length(tips))
tips_n[i] = which(phy$tip.label == tips[i])
rows = vector("numeric",length(tips_n))
for (i in 1:length(tips_n))
rows[i] = which(phy$edge[,2] == tips_n[i])
nodes = sort(phy$edge[rows[which.min(rows)]:rows[which.max(rows)],1])
return(nodes[1])
} else if (is.numeric(tips) == TRUE){
rows = vector("numeric",length(tips))
for (i in 1:length(tips))
rows[i] = which(phy$edge[,2] == tips[i])
nodes = sort(phy$edge[rows[which.min(rows)]:rows[which.max(rows)],1])
return(nodes[1])
} else if (is.null(tips) == TRUE){
stop('tips vector is empty')
}
}
Liam, I actually wrote this to "improve" your collapse.to.star script in
phytools. I'm posting it below.
It basically allows collapsing N number of nodes at once given a tree and a
vector of nodes (node) or a single node.
Could you give it a try and let me know if its buggy?
I worked fine for me, with the new which.node function.
Cheers, and thanks again!
Santiago
collapse.to.star2 <- function (tree, node)
{
require(phangorn)
require(phytools)
which.node <- function(phy, tips)
{
if (is.character(tips) == TRUE){
tips_n = vector("numeric",length(tips))
for (i in 1:length(tips))
tips_n[i] = which(phy$tip.label == tips[i])
rows = vector("numeric",length(tips_n))
for (i in 1:length(tips_n))
rows[i] = which(phy$edge[,2] == tips_n[i])
nodes = sort(phy$edge[rows[which.min(rows)]:rows[which.max(rows)],1])
return(nodes[1])
} else if (is.numeric(tips) == TRUE){
rows = vector("numeric",length(tips))
for (i in 1:length(tips))
rows[i] = which(phy$edge[,2] == tips[i])
nodes = sort(phy$edge[rows[which.min(rows)]:rows[which.max(rows)],1])
return(nodes[1])
} else if (is.null(tips) == TRUE){
stop('tips vector is empty')
}
}
o_tree <- tree
if (is.vector(node) == FALSE){
if (is.numeric(node) == TRUE){
tt <- splitTree(tree, split = list(node = node, bp =
tree$edge.length[which(tree$edge[,2] == node)]))
ss <- starTree(species = tt[[2]]$tip.label, branch.lengths =
diag(vcv(tt[[2]])))
ss$root.edge <- 0
tree <- paste.tree(tt[[1]], ss)
return(tree)
} else {
stop('node is not numeric')
}
} else if (is.vector(node) == TRUE){
if (is.numeric(node) == TRUE){
subtree_list = vector("list", length(node))
for (i in 1:length(node)){
if (is.null(subtree_list[[1]]) == TRUE){
tt <- splitTree(tree, split = list(node = node[1], bp =
tree$edge.length[which(tree$edge[,2] == node[1])]))
ss <- starTree(species = tt[[2]]$tip.label, branch.lengths =
diag(vcv(tt[[2]])))
ss$root.edge <- 0
subtree_list[[1]] <- ss
tree <- paste.tree(tt[[1]], ss)
} else {
new_node <- which.node(tree, o_tree$tip.label[Descendants(o_tree,
node[i])[[1]]])
tt <- splitTree(tree, split = list(node = new_node, bp =
tree$edge.length[which(tree$edge[,2] == new_node)]))
ss <- starTree(species = tt[[2]]$tip.label, branch.lengths =
diag(vcv(tt[[2]])))
ss$root.edge <- 0
#subtree_list[[i]] <- ss
tree <- paste.tree(tt[[1]], ss)
}
}
return(tree)
} else {
stop('node vector is not numeric')
}
}
}
2013/11/25 Liam J. Revell >
> Hi Santiago.
>
> You should check out findMRCA in the phytools package. It does what I
> think you want to do. Also, I encountered cases in which which.node doesn't
> seem to work properly. One difference is that which.node can
> (theoretically) take tip labels or numbers as input. To use tip numbers as
> input in findMRCA you could just do:
>
> x<-findMRCA(tree,tree$tip.label[tip.numbers])
>
> since tip numbers are just the order of tree$tip.label.
>
> I believe that findMRCA does not have any serious bugs. Please let me know
> if you find that this is not the case. It uses (fastMRCA which uses)
> Ancestors from phangorn internally.
>
> All the best, Liam
>
> Liam J. Revell, Assistant Professor of Biology
> University of Massachusetts Boston
> web: http://faculty.umb.edu/liam.revell/
> email: liam.rev...@umb.edu 'liam.rev...@umb.edu');>
> blog: http://blog.phytools.org
>
>
> On 11/25/2013 4:10 PM, Santiago Sánchez wrote:
>
>> which.node <- function(phy, tips)
>> {
>> if (is.character(tips) == TRUE){
>> tips_n = vector("numeric",length(tips))
>> for (i in 1:length(tips))
>> tips_n[i] = which(phy$tip.label == tips[i])
>> rows = vector("numeric",length(tips_n))
>> for (i in 1:length(tips_n))
>
[R-sig-phylo] which.node function
Hello list,
I'm not aware that a function such as which.node exists in any phylo
package available. Or at least one that can retrieve the node number from a
tree (phy) given a list of tips (tips). For this reason I went ahead and
wrote a simple function (see below). The function allows tips to be given
in "character" or "numeric" classes.
If someone finds it useful or so if any of the developers would like to
make it part of their packages, go ahead.
which.node <- function(phy, tips)
{
if (is.character(tips) == TRUE){
tips_n = vector("numeric",length(tips))
for (i in 1:length(tips))
tips_n[i] = which(phy$tip.label == tips[i])
rows = vector("numeric",length(tips_n))
for (i in 1:length(tips_n))
rows[i] = which(phy$edge[,2] == tips_n[i])
nodes = sort(phy$edge[rows,1])
return(nodes[1])
} else if (is.numeric(tips) == TRUE){
rows = vector("numeric",length(tips))
for (i in 1:length(tips))
rows[i] = which(phy$edge[,2] == tips[i])
nodes = sort(phy$edge[rows,1])
return(nodes[1])
} else if (is.null(tips) == TRUE){
stop('tips vector is empty')
}
}
All the best,
Santiago
--
Santiago Sánchez-Ramírez
Department of Ecology and Evolutionary Biology, University of Toronto
Department of Natural History (Mycology), Royal Ontario Museum
100 Queen's Park
Toronto, ON
M5S 2C6
Canada
Other email: santiago.sanc...@mail.utoronto.ca
Tel. 416-586-8025
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
