This is cool — so this looks like it will work on a set of posterior trees even if they don’t all share the same topology? Jake
> On Apr 10, 2018, at 3:37 AM, jschenk <jsch...@georgiasouthern.edu> wrote: > > Thank you Graham, Liam, and Emmanuel for your suggestions. Applying the > .uncompressTipLabel did the trick for my code. If anyone is interested in > using the code to identify the 95% HPD for a single node from a set of > chronograms, you can find the code here > (https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R > <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R><https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R > <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R>>) or > below. > > ~John > > > > library(ape) > library(coda) > > #Function to estimate the age of a node given two members of a clade. > AgeDensity <- function(phy, species1, species2){ > NodeNumber <- mrca(phy)[species1, species2] > ages <- branching.times(phy)[as.character(NodeNumber)] > return(as.numeric(ages)) > } > > #read in tree > tree <- read.nexus(file="") > > #Apply function across the posterior distribution to obtain node ages. This > will take some time to run, depending on the number of trees. It took me > about 20 minutes. > NodeAges <- unlist(lapply(.uncompressTipLabel(tree), AgeDensity, "species1", > "species2")) > > #estimate the HPD interval for the node ages > HPD <- HPDinterval(as.mcmc(NodeAges), prob = 0.95) > > #Density plot for 95% HPD > plot(density(as.numeric(NodeAges[which(NodeAges > HPD[1, 1] & NodeAges < > HPD[1, 2])])), main="", las=1, xlab = "Million of years before present", lwd > = 2) > > > ______________________________________________________________ > John J. Schenk, Ph.D. > Assistant Professor of Plant Biology > Georgia Southern University Herbarium (GAS), Curator > Department of Biology > 4324 Old Register Road > Georgia Southern University > Statesboro, GA 30460-8042 > Office: 2260 Biology Building > Office phone: (912) 478-0848 > Lab website: sites.google.com/a/georgiasouthern.edu/schenk > <http://sites.google.com/a/georgiasouthern.edu/schenk> > Herbarium website: sites.google.com/a/georgiasouthern.edu/gasherbarium > <http://sites.google.com/a/georgiasouthern.edu/gasherbarium> > > > > > > > >> On Apr 9, 2018, at 10:56 AM, Liam J. Revell <liam.rev...@umb.edu> wrote: >> >> Dear John. >> >> You could try running .uncompressTipLabel on the "multiPhylo" object. Let us >> know if that works. >> >> All the best, Liam >> >> Liam J. Revell, Associate Professor of Biology >> University of Massachusetts Boston >> & Profesor Asociado, Programa de Biología >> Universidad del Rosario >> web: http://faculty.umb.edu/liam.revell/ >> >> On 4/9/2018 9:46 AM, jschenk wrote: >>> Hi Folks, >>> I have been banging my head against what appears to be an easy coding >>> problem for a while now and haven’t been able to hack my way out of it. I >>> am running a function to identify a posterior set of node ages for a >>> particular node. The function I wrote works just fine, but when I use >>> lapply to sample across a posterior distribution, I get the “not of class >>> ‘phylo’” error, even when I input a single tree of class phylo. I have >>> tried every typical solution (e.g., reassigning class), but haven’t >>> identified a solution. Does anyone know a workaround? >>> Please see the example code below. >>> Thanks, >>> John >>> library(ape) >>> #simulate a single tree with 20 tips >>> simtree <- rtree(20) >>> #Make sure the tree exists >>> plot(simtree) >>> #Function I wrote to find a node and then tell me the age of the node. I >>> realize that the simulated tree is not ultrametric in this example, in real >>> life it will be - ultrametric trees also result in the same error >>> AgeDensity <- function(phy, species1, species2){ >>> NodeNumber <- mrca(phy)[species1, species2] >>> ages <- branching.times(phy)[as.character(NodeNumber)] >>> return(as.numeric(ages)) >>> } >>> #check the class of the tree object, it will say that it is of class phylo >>> class(simtree) >>> #Run my function AgeDensity and it works just fine >>> AgeDensity(simtree, "t3", "t15") >>> #When I use the lapply function, I get an error that the object is not a of >>> class phylo, although I already verified that it is of class phylo. >>> lapply(simtree, AgeDensity, species1="t3", species2="t15") >>> #here is the same analysis conducted with multiple trees >>> multiTrees <- rmtree(20, 10) >>> class(multiTrees) >>> #I get the same error when I run my function across multiple trees >>> lapply(multiTrees, AgeDensity, species1="t3", species2="t15") >>> ______________________________________________________________ >>> John J. Schenk, Ph.D. >>> Assistant Professor of Plant Biology >>> Georgia Southern University Herbarium (GAS), Curator >>> Department of Biology >>> 4324 Old Register Road >>> Georgia Southern University >>> Statesboro, GA 30460-8042 >>> Office: 2260 Biology Building >>> Office phone: (912) 478-0848 >>> Lab website: sites.google.com/a/georgiasouthern.edu/schenk >>> Herbarium website: sites.google.com/a/georgiasouthern.edu/gasherbarium >>> [[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/ > > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-phylo mailing list - R-sig-phylo@r-project.org > <mailto:R-sig-phylo@r-project.org> > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> > Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ > <http://www.mail-archive.com/r-sig-phylo@r-project.org/> [[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/
