[R-sig-phylo] Problem with make.simmap in Phytools

Liam J. Revell

2018-11-01 03:49:51 UTC

Dear Kimberly.

I see one error in your script which is that when you populate the list
of sets of stochastic map trees you call the object 'simtrees' - but
when you try to summarize the results, you attempt to pass an object
called 'simtree' (no 's') to the summary method.

A couple of other comments.

(1) The function call list(500) creates a list with one object, the
numeric value 500, rather than a list with 500 elements.

(2) In general, it is not necessary to iterate tree by tree over you
"multiPhylo" object in the way you have. Instead, the simple call:

simtrees<-make.simmap(trees,trait,model,nsim)

will produce a "multiSimmap" object containing nsim x length(trees)
stochastic map trees.

(3) Finally, you can call the S3 summary method for the object class
"multiSimmap" directly on this single, "multiSimmap" object containing
all the stochastic map trees - rather than on each set of nsim trees in
turn. If you decide to do that, however, you should keep in mind that
trees from a Bayesian posterior sample may differ one from the other in
topology. For this I recommend using the argument check.equal OR
supplying a reference tree via the argument ref.tree. In the former
case, the method first checks if all trees are equal, and if they are
not, computes a consensus tree and the posterior probabilities at each
node of the consensus topology conditional on that node being present in
each sampled tree. If ref.tree is supplied, the method will do the same
but using the reference tree instead of a consensus tree. In the case of
Bayesian MCMC, the reference tree could be (for instance) the MCC tree.
Useful examples may be found be searching for "ref.tree" on my blog.

I hope that this helps you find your error. If it does not, please free
to send me your saved workspace & I will try to figure it out.

All the best, Liam

Liam J. Revell
Associate Professor, University of Massachusetts Boston
Profesor Asistente, Universidad Católica de la Ssma Concepción
web: http://faculty.umb.edu/liam.revell/, http://www.phytools.org

Post by Kimberly Lynn Foster
Hello, I am currently trying to run make.simmap in Phytools ver 0.6-60 and Geiger ver 2.0.6, I am running the function on a subsample of 500 trees from the posterior distribution. I have included the full code, but there error I get is found in the last line of code
simtrees<-list(500)
YY<-list(500)
for (i in 1:500){
simtrees[[i]]<-make.simmap(CrenuchidTreesLadderized[[i]], habitattrait, model="ER", nsim=1000)
YY[[i]]<- describe.simmap(simtree[[i]], plot=F)
}
no applicable method for 'isSymmetric' applied to an object of class "c('double', 'numeric')"
require(phytools)
require(geiger)
habitat<-read.csv("C:/Users/XXX/Dropbox/Crenuchid R Files/Crenuchid R Files/Crenuchid_2state_BvsP.csv", row.names=1)
habitat[habitat==0] <- "Pelagic"
habitat[habitat==1] <- "Benthic"
habitat_color <- habitat
habitat_color[habitat=="Benthic"] <- "#6EA0FF" ### blue
habitat_color[habitat=="Pelagic"] <- "#FFB322" ### Orange
############################# simmap multiple trees ################################################.
alltrees<-read.nexus("posteriorsamples.tre.trprobs") #reading trees from posteior distribution
wib <- name.check(phy=alltrees[[1]], data=habitat)
wib ## will show what is in tree not data, data not tree
# ??name.check
# Now to drop tips
length(alltrees)
prunedpostsamp<-list(500)
for (i in 1:500){
prunedpostsamp[[i]]<-drop.tip(alltrees[[i]], wib$tree_not_data)
}
wib <- name.check(phy=prunedpostsamp[[1]], data=habitat)
prunedpostsamp[[1]]$tip.label
#Let's ladderize all trees in the distribution with lapply() and plot the first tree again
CrenuchidTreesLadderized <- lapply(prunedpostsamp, ladderize)
plot(CrenuchidTreesLadderized[[1]], cex=0.4) #Plot the first ladderized tree
#However, does the order of tips in the plot match the order of tiplabels in the tree object? #### NO help? is this the issue ???
#Please compare######################################################################################
CrenuchidTreesLadderized[[1]]$tip.label
habitattrait<-habitat[CrenuchidTreesLadderized[[1]]$tip.label,]
names(habitattrait)<- CrenuchidTreesLadderized[[1]]$tip.label
length(habitattrait) ## 54
class(habitattrait) ## character
#Compare tree with data
compare <- treedata(CrenuchidTreesLadderized[[1]], habitattrait, sort=TRUE) #R will inform you of any inconsistencies and how they were solved!
?treedata
####### SIMMAP for multiple trees this can take a while while while ########
simtrees<-list(500)
YY<-list(500)
for (i in 1:500){
simtrees[[i]]<-make.simmap(CrenuchidTreesLadderized[[i]], habitattrait, model="ER", nsim=1000)
YY[[i]]<- describe.simmap(simtree[[i]], plot=F)
}
[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/